Sensitivity analysis and constraint of the parameter space of Earth System configuration of JULES

1 Introduction

Update of explore-JULES-ES-1p0.Rmd that uses the new packages and functions

We build and test gaussian process emulators, perform a sensitivity analyis, and constrain the input space of Jules.

Andy thinks there might be mileage in analysing the atmospheric growth, which is here. /home/h01/hadaw/Research/ES_PPE/Oct20

At the moment, this vignette is hampered by the fact that emulators are failing on a few of the outputs which represent change over the historical period. The emulator is fine for predicting absolute values in the modern period.

Andy would like to see timeseries of: cVeg, cSoil and nbp, npp in GtC and GtC/yr.

1.1 Preliminaries

Load libraries, functions and data.

1.2 Carbon budget data

Section 2.5 in Friedlingstein et al. describes how the land carbon sink is estimated.

1.2.1 Variables available for analysis

Load a single ensemble member to examine contents

1.3 Plot timeseries of aggregated global timeseries

1.4 Where did the model fail to run?

(“probability of run failure” Could be an interesting project - logit transformation for probability of run failure? Some other ML like SVM?

There are clear run failure thresholds in the parameters rootd_ft_io and lai_max, and quite strong visual indications that a_wl_io and bio_hum_cn matter.

NEXT, histograms of failure with the LHS value

Could try an SVM, logistic regression or similar.

1.5 Testing Gaussian Process emulators

1.5.1 NPP as a test

First, use mean NPP as an example. How does NPP respond to each parameter? NAs are removed, but zero values are still included.

1.6 Relationship between modern NPP and change since 1850.

Some outputs (e.g fLuc, fHarvest) have an almost perfect 1:1 relationship between modern value and change, some (nbp, npp, treeFrac) quite or moderately strong, and some (csoil, cveg) very weak or non-existant.

1.7 A clear threshold in the F0 parameter for NPP

It appears that this ensemble is less “clear cut” in having an output that clearly distinguishes between “failed” (or close to it), and “not failed”.

Having said that, having an F0 over a threshold seems to kill the carbon cycle, as before. Here, we’ve set a threshold of 0.9 (on the normalised scale) for F0, and we remove members of the ensemble with a larger F0 than that when we build emulators.

[1] \nbp_lnd_sum\
[1] \Carbon Mass Flux out of Atmosphere Due to Net Biospheric Production on Land\
[1] \year\
[1] \year\
[1] \nbp_lnd_mean\
[1] \Carbon Mass Flux out of Atmosphere Due to Net Biospheric Production on Land\
[1] \fLuc_lnd_sum\
[1] \Net Carbon Mass Flux into Atmosphere Due to Land-Use Change\
[1] \fLuc_lnd_mean\
[1] \Net Carbon Mass Flux into Atmosphere Due to Land-Use Change\
[1] \npp_nlim_lnd_sum\
[1] \Net Primary Production on Land as Carbon Mass Flux\
[1] \npp_nlim_lnd_mean\
[1] \Net Primary Production on Land as Carbon Mass Flux\
[1] \cSoil_lnd_sum\
[1] \carbon mass in model soil pool\
[1] \cSoil_lnd_mean\
[1] \carbon mass in model soil pool\
[1] \cVeg_lnd_sum\
[1] \carbon mass in vegetation\
[1] \cVeg_lnd_mean\
[1] \carbon mass in vegetation\
[1] \landCoverFrac_lnd_sum\
[1] \Plant Functional Type Grid Fraction\
[1] \landCoverFrac_lnd_mean\
[1] \Plant Functional Type Grid Fraction\
[1] \fHarvest_lnd_sum\
[1] \Carbon Mass Flux into Atmosphere Due to Crop Harvesting\
[1] \fHarvest_lnd_mean\
[1] \Carbon Mass Flux into Atmosphere Due to Crop Harvesting\
[1] \lai_lnd_sum\
[1] \leaf area index\
[1] \lai_lnd_mean\
[1] \leaf area index\
[1] \rh_lnd_sum\
[1] \Total Heterotrophic Respiration on Land as Carbon Mass Flux\
[1] \rh_lnd_mean\
[1] \Total Heterotrophic Respiration on Land as Carbon Mass Flux\
[1] \treeFrac_lnd_sum\
[1] \Fractional cover of each surface type\
[1] \treeFrac_lnd_mean\
[1] \Fractional cover of each surface type\
[1] \c3PftFrac_lnd_sum\
[1] \Percentage Cover by C3 Plant Functional Type\
[1] \c3PftFrac_lnd_mean\
[1] \Percentage Cover by C3 Plant Functional Type\
[1] \c4PftFrac_lnd_sum\
[1] \Percentage Cover by C3 Plant Functional Type\
[1] \c4PftFrac_lnd_mean\
[1] \Percentage Cover by C3 Plant Functional Type\
[1] \shrubFrac_lnd_sum\
[1] \Percentage Cover by Shrub\
[1] \shrubFrac_lnd_mean\
[1] \Percentage Cover by Shrub\
[1] \baresoilFrac_lnd_sum\
[1] \Bare Soil Percentage Area Coverage\
[1] \baresoilFrac_lnd_mean\
[1] \Bare Soil Percentage Area Coverage\
[1] \residualFrac_lnd_sum\
[1] \Percentage of Grid Cell That Is Land but neither Vegetation Covered nor Bare Soil\
[1] \residualFrac_lnd_mean\
[1] \Percentage of Grid Cell That Is Land but neither Vegetation Covered nor Bare Soil\

1.8 Re-examine bwl_io

A closer look at bwl_io now that the impact of f0_io has been removed shows a large number of “zero” NPP when bwl_io is at low values, which could well be the source of apparent sensitivity.

Take only higher values of bwl_io for the next round of constraints.

2 History matching and generation of new ensemble members

To use History Matching, we need to specify targets for various model outputs. We treat these targets as “observations” with an uncertainty where a model run marked as “implausible” (beyond 3sd) matches the hard boundaries previously identified by A. Wiltshire as being desirable.

Choose the centre of the (implied) uniform distribution. cs_gb.target = (3000 - 750) / 2 = 1125 cv.target = (800 - 300) / 2 = 250 npp_n_gb.target = (80 - 35) / 2 = 22.5

nbp.target = 0 (gpp.target = 75) (runoff.target = 1)

(to do: visualise implausibility of design and loo emulated values of design as histogram)

2.1 visualise the constraints

2.2 Timeseries that match the constraints.

Plot the original ensemble, plus the stuff that matches the constraints. Include anomalies. “Level 2” constraints are broad constraints produced by AW

2.3 Pairs plots of level 2 input space

These are inputs that match the level 2 constraints. Only 37 of the original 500 members match these constraints. It’s noticable how wide a region of parameter space these still come from.

We find that adding a constraint to level 2 constrains the behaviour of the soil carbon pool, but not the vegetation carbon pool. Probably not a surprise, given the Soil Carbon pool (and changes) are much larger than the vegetation carbon pool.

2.4 Emulated density pairs plot of input space that matches level 2 constraints

2.5 Constrain with Cumulative NBP

Email from ER:

Hello,

I’ve had a look for some multi-model ranges…

From the ILAMB CMIP6 pages. NBP accumulated between 1959 and 2014: (227 INM) -16.3 Gt to +106 Gt GPP 1980-2014 mean: 104-144 Gt/yr (157 INM, 159 BCC) cVeg 1996-2006 mean: 338-582 Gt (635 INM) cSoil 2000-2001 mean: 319-1760 Gt

From Spencer’s recent paper about a subset of CMIP6 models. NBP 2000-2009 mean: 0.45-2.24 Gt/yr

From Trendy 2019 ensemble, all are 1985-2014 means. NBP: 0.21-1.82 Gt/yr GPP: 107-165 Gt/yr cVeg: 335-869 Gt cSoil: 626-3595 Gt

I can update the Trendy results to the latest 2020 version and/or change the meaning period if needed.

Cheers, Eddy.

If we constrain on cumulative NBP that matched ILAMB 1959 - 2013 (the latest year in the record, we choose those ensemble members which have the lowest cumulative NBP. Many ensemble members have a higher cumulative NBP.)

newdata not named: newdata's variables are inherited from the designnewdata not named: newdata's variables are inherited from the designnewdata not named: newdata's variables are inherited from the designnewdata not named: newdata's variables are inherited from the design
[1] 13.25

2.6 How do we find outputs that are inconsistent?

It’s clear that something rules out all the “good” values of cumulative_nbp. Which output is it?

1. There is a very strong relationship (almost 1:1) between nbp_lnd_sum and cumulative_nbp. The cumulative NBP constraint keeps the ‘modern’ nbp low, while Andy’s constraints do not.

2. All the cumulative NBP-constrained ensemble members have low modern vegetation carbon

3. Modern NPP and soil carbon #

2.7 Augment the design.

The function addNroyDesignPoints builds an emulator for each model output in Y. It compares the output of each emulator at a number of candidate desin points, and chooses a space-filling set of them that that are Not Ruled Out Yet (statistically close to the observation at Y_target).

2.8 Write the augmented design

The function write_jules_design here simply takes the points calculated by addNroyDesignPoints and writes them to configuration files.

2.9 Check the design

Check that the augmented design produces what we expect. New ensemble members should be somewhat constrained within the boundaries of the original design, if the comparison to data offers any constraint.

Visualise the predicted outputs at the NROY points of the old design, and the new suggested design.

2.10 Visualising the emulated outputs at the proposed new design points.

Interestingly, these aren’t perfectly within the original hard boundaries set by Andy, even though I’ve set those boundaries to be the +- 4 standard deviation threholds in the History Match. I suggest this is because there is model discrepancy, and that there is considerable wriggle room induced from emulator uncertainty.

In particular, it appears that vegetation carbon is difficult to keep high, and that many NROY proposed members have a fairly low vegetation carbon. This might need a discrepancy term, or adjusting in some other way. It certainly needs exploring, and a OAAT plot might give clues as to the parameters to choose.

2.11 Atmospheric growth

We can calculate the atmospheric growth as Fossil fuel emissions (GCB) - Ocean uptake (GCB) - net land flux (JULES NBP)

Conversely, we can calculate the land uptake as NBP = FF - ocean - ag (Do we have AG in JULES?)

some values will be clipped
null device 
          1 

2.11.1 another way to calclulate atmospheric growth

These two methods give very similar outcomes, bar an outlier.

2.11.2 Land carbon Sink?

I think to make a fair comparison, we might need to take the diff of the anomaly (or just the timeseries), as the above represents the total change since 1850, and I think all the “land sink” stuff looks at the year-to-year changes.

2.11.3 Identify parts of parameter space where the land carbon sink is at its maximum

What do these parts of parameter space tell us about what we need to do to have the model reproduce reality?

2.11.4 what parts of parameter space are the closest?

How good are the four emulators that we’ve built? Are there biases? (there’s no real evidence of this)

What do the emulators make of the design points which actually make Andy’s “hard boundary” criteria? If we leave them out, do they still place the output within the hard boundaries?

2.12 Check the emulators that produce the new design

Do a leave-one-out cross validation of points inside the hard boundaries using the wave1 fits.

Error: object 'pred_unif' not found

We see in the leave-one-out analysis that the emulator is consistently under-predicting the vegetation carbon (though the uncertainty estimate often covers the actual value).This suggests (1) that there isn’t really a huge problem with a model discrepancy (or at least that isn’t the only problem), and (2) the history matching is working as it should, and taking into account a not-great emulator.

2.13 Compare the straight km with a two-step emulator

Is a two-step emulator any better at emulating those crucial points which fall within aw’s hard boundaries? First, create a list of emulator fits.

It doesn’t appear that the two-step emulator (here plotted in red) is doing any better than the regular emulator.

Two other things I can think of to check: 1) how about using “multistart” to choose different starting conditions for optimising the emulators and 2) Using a flat prior for the mean function rather than a linear prior.

3 Exploring new constraints

Input from Eddy Robertson:

“We can assume that the majority of vegetation carbon is stored in tree trunks so, the carbon density of trees is approximately (cVeg_lnd_mean / treeFrac_lnd_mean). Although I suppose it’s possible that in the PPE the shrubs have become tree-like, so it might be interesting to plot cVeg_lnd_mean versus (shrubFrac_lnd_mean + treeFrac_lnd_mean) as well. The first question is whether cVeg_lnd_mean increases linearly with treeFrac_lnd_mean. Then if this turned out to be interesting, we could produce a treeCVeg_lnd_mean output.”

Is the range of carbon densities narrower in the “hard boundary” set?

3.1 Now that we have constraints, we can explore their effect on the output

3.2 How much does constraining an output constrain other outputs?

For example, what effect has the AW “hard boundary” constraint produced on ranges of the other (marginal) outputs? Which constraint is doing the majority of the work?

A function to find the proportion of output space that is removed (or retained) when applying a constraint.

All constraints together take the normalized range of all the outputs to a little under a half of the original range, on average. Of course, four of the inputs are constrained here, which has quite a large impact.

We can see what impacts the constraints have individually

The NBP constraint has no marginal impact

Vegetation carbon (cVeg) has the strongest impact on everything else. This is no suprise, given that so few of the ensemble members seem to simulate carbon vegetation high enough.

3.3 Writing robust, parallelised, multistart emulators

4 Direct twoStep emulator using foreach