GCEL / CARDAMOM

The CARbon DAta MOdel FraMework. Computer software that retrieves terrestrial carbon (C) cycle variables by combining C cycle observations with a mass balance model.
GNU General Public License v3.0
2 stars 0 forks source link

Setting new parameter priors (for leaf lifespan) #28

Open DTMilodowski opened 1 year ago

DTMilodowski commented 1 year ago

Hi all,

I'm going to try setting priors for the leaf lifespan for the forest classes in the DARE-UK analysis, but have a couple of questions.

The first is whether it is sufficient to set the prior directly on the lifespan parameter (for DALEC-CDEA-ACM2)?

The second is which of the R functions needs to be edited to add in this information into the parameter prior list?

My initial implementation will be on the old version of CARDAMOM (for legacy reasons) but I'll eventually implement in the main dev branch too. I will also keep notes to document the process for similar changes in the future.

Thanks for the feedback!

DTMilodowski commented 1 year ago

ok, some quasi-intelligent searching reveals that parameter priors can be set here in the binary_data.r function

Which just leaves the first question - whether leaf lifespan can be set directly on the parameter?

Thanks!

lsmallma commented 1 year ago

Well searched!

I’d suggest setting it as an other prior which can then up used in the model_likelihood.f90 file.

I can check which have not been used before tomorrow morning.

Thanks, Luke

Dr T. Luke Smallman Research Fellow (NCEO)

School of GeoSciences Crew Building University of Edinburgh EH9 3JN, UK Edinburgh

On 29 Nov 2022, at 15:41, David Milodowski @.***> wrote:

This email was sent to you by someone outside the University. You should only click on links or attachments if you are certain that the email is genuine and the content is safe.

ok, some quasi-intelligent searching reveals that parameter priors can be set herehttps://github.com/GCEL/CARDAMOM/blob/00815799bfad9738930196c4b8e3ec5dd38be011/R_functions/binary_data.r#L226 in the binary_data.rhttps://github.com/GCEL/CARDAMOM/blob/00815799bfad9738930196c4b8e3ec5dd38be011/R_functions/binary_data.r function

Which just leaves the first question - whether leaf lifespan can be set directly on the parameter?

Thanks!

— Reply to this email directly, view it on GitHubhttps://github.com/GCEL/CARDAMOM/issues/28#issuecomment-1330844902, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AGDTODIWPK64LLM5QR7SIEDWKYPZTANCNFSM6AAAAAASOTXWK4. You are receiving this because you are subscribed to this thread.Message ID: @.***>

The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336. Is e buidheann carthannais a th’ ann an Oilthigh Dhùn Èideann, clàraichte an Alba, àireamh clàraidh SC005336.

DTMilodowski commented 1 year ago

Just adding here for future reference, the recommendation to set the leaf life span prior constraint as an "other prior" rather than directly on the "leaf lifespan" parameter (P5 in DALEC-CDEA-ACM2) is that this parameter sets the baseline turnover fraction in each timestep, but the actual effective leaf lifespan in the model depends on other components too. Thus, the better application for the desired prior constraint is through calculation of the mean transit time in the calibrated model vs. the prior leaf lifespan estimate. Thanks @lsmallma for the clarification!

There are additional potential complexities around how this (and other priors) might be weighted relative to individual observations, but that is part of a broader set of challenges in how to appropriately weight the uncertainties from diverse observation streams.

DTMilodowski commented 1 year ago

Hi all. I'm still testing, but thought I'd post my initial approach here for future reference (@tim-green15)

First, I added some clauses into my control file. In the future, we could have specific scripts that are called internally within the R framework.

## Set some prior parameters...
set_leaflifespan_prior = FALSE
if (lc_class == 'coniferous_forest'){
    set_leaflifespan_prior = TRUE
    leaflifespan_prior = 3.05
    leaflifespan_prior_unc = 0.81
}
if (lc_class == 'broadleaf_forest'){
    set_leaflifespan_prior = TRUE
    leaflifespan_prior = 0.53
    leaflifespan_prior_unc = 0.15
}

Then, I added this snippet into the extract_obs.r script:

    ###
    ## Get leaf lifespan information
    if (set_leaflifespan_prior == TRUE) {
        print(paste("\tsetting leaf lifespan prior = ",leaflifespan_prior,"+/-",leaflifespan_prior_unc,sep=" "))
        leaflifespan=leaflifespan_prior
        leaflifespan_unc=leaflifespan_prior_unc
    } else {
        # assume no data available
        leaflifespan=-9999 ; leaflifespan_unc=-9999
    }

I also added this line into the binary_data.r script:

      OTHERPRIORS[7] = OBS$leaflifespan ; OTHERPRIORUNC[7] = OBS$leaflifespan_unc # leaf lifespan

Note that slot 7 is currently open across all the UoE models, and should now be considered as allocated to the "Leaf lifespan" slot in future development.

And then in the fortran code, I updated the MODEL_LIKELIHOOD.f90 file to add the phenological leaf lifespan (excluding disturbance fluxes) vs. prior and uncertainty. Obviously, any added variables are predeclared earlier in the function.

    ! Check leaf lifespan.
    ! NOTE: this arrangement explicitly neglects the impact of disturbance on
    ! residence time (i.e. no fire and biomass removal)
    if (DATAin%otherpriors(7) > -9998) then
        ! Estimate the mean annual input to the wood pool (gC.m-2.day-1) and
        ! remove the day-1 by multiplying by residence time (day)
        ! Foliage MRT (years-1)
        tot_exp = 0
        fol_hak = 0
        do n = 1, DATAin%nodays
            if (DATAin%M_POOLS(n,2)>0) then
               fol_hak = fol_hak + 1
               tot_exp = tot_exp + (DATAin%M_FLUXES(n,10) / DATAin%M_POOLS(n,2))
            end if
        end do
        tot_exp = tot_exp / fol_hak
        tot_exp = 1 / (tot_exp * 365.25)
        likelihood = likelihood - ((tot_exp - DATAin%otherpriors(7)) / DATAin%otherpriorunc(7))**2
    endif

Hopefully that will do the trick. Still a work in progress and currently being tested.