You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Calibration of a new community type within DVM-DOS-TEM often involves tuning vegetation and soil carbon related parameters to match observations or literature compilations of carbon stocks and fluxes. While such calibration is necessary to begin representing ecosystem dynamics, it is still only part of the picture. In DVM-DOS-TEM, the thermal and hydrological properties of the soil exert considerable influence on carbon stocks and fluxes through processes affected by soil moisture and temperature. It is recommended that the thermo-hydrological parameter space is explored during calibration, to best represent important processes affecting ecosystem dynamics.
This example provides a framework for running sensitivity analysis on parameters with the greatest influence on the thermal and hydrological behavior of the model. In this case, we will be comparing model output soil moisture and temperature with site-based observations, in an attempt to determine an optimal set of parameters.
A list of important parameters with units, ranges, and literature references was developed by Helene Genet.
Parameters influencing soil hydro/thermal state
cmt_envground.txt
tcsolid(m): Soil thermal conductivity for moss (W/mK)
current value: 0.25, range: 0.005-0.5
Ekici et al. 2015,Jiang et al. 2015, O’Donnell et al. 2009
tcsolid(f): Soil thermal conductivity for fibric (W/mK)
current value: 0.25, range: 0.005-0.5
Ekici et al. 2015,Jiang et al. 2015, O’Donnell et al. 2009
tcsolid(h): Soil thermal conductivity for humic (W/mK)
current value: 0.25, range: 0.02-2.0
Ekici et al. 2015,Jiang et al. 2015, O’Donnell et al. 2009
porosity(m): porosity for moss layers (m3/m3)
current value: 0.98, range: 0.85-0.99
O’Donnell et al. 2009
porosity(f): porosity for fibric layers (m3/m3)
current value: 0.95, range: 0.85-0.99
O’Donnell et al. 2009
porosity(h): porosity for humic layers (m3/m3)
current value: 0.8, range: 0.7-0.9
O’Donnell et al. 2009
bulkden(m): bulk density for moss (g/m3)
current value: 25,000, range: 10,000 - 80,000
Tuomi et al. 2020, Rodionov et al. 2007, O’Donnell et al. 2009
bulkden(f): bulk density for fibric (g/m3)
current value: 51,000, range: 20,000 - 200,000
Tuomi et al. 2020, Rodionov et al. 2007, O’Donnell et al. 2009
bulkden(h): bulk density for humic (g/m3)
current value: 176,000, range: 100,000 - 800,000
Tuomi et al. 2020, Rodionov et al. 2007, O’Donnell et al. 2009
hksat(m): hydraulic conductivity at saturation for moss (mm/s)
current value: 0.15, range: 0.0002 - 30
Ekici et al. 2015, Letts et al. 2000, Liu et al. 2019
hksat(f): hydraulic conductivity at saturation for fibric (mm/s)
current value: 0.28, range: 0.0002 - 30
Ekici et al. 2015, Letts et al. 2000, Liu et al. 2019
hksat(h): hydraulic conductivity at saturation for humic (mm/s)
current value: 0.002, range: 0.00004 - 2.01
Ekici et al. 2015, Letts et al. 2000, Liu et al. 2019
nfactor(s): Summer nfactor
current value: 1.5, range: 0.2-2.0
Kade et al. 2006, Klene et al. 2001
snwalbmax
current value: 0.8, range: 0.7-0.85
Te Beest et al. 2016, Petzold et al. 1975, Loranty et al. 2011
snwalbmin
current value: 0.4, range: 0.4-0.6
Te Beest et al. 2016, Petzold et al. 1975, Loranty et al. 2011
cmt_dimground.txt
snwdenmax (kg/m3)
current value: 250, range: 100-800
Domine et al. 2016, Gerland et al. 1999, Muskett 2012
snwdennew (kg/m3)
current value: 50, range: 10-250
Domine et al. 2016, Gerland et al. 1999, Muskett 2012
cmt_bgcsoil.txt
rhq10
current value: 2, range: 1.6-2.4
rhq10_w
current value: 2, range: 0.85-0.99
SA setup
We will use a few useful imports from the 'scripts' directory in addition to a few other packkages.
import sys
import os
# setting path
sys.path.append('/work/scripts')
import pandas as pd
import seaborn as sns
import Sensitivity as sa
import output_utils as ou
Define working directory and run setup string for the sensitivity driver.
The sensitivity driver has the option to define a percent to perturb the parameters, however this does not work so well for parameters spanning multiple orders of magnitude, so we will manually specify the bounds. Additionally, we will specify which parameters should be sampled on a log-uniform scale; i.e. the parameters spanning multiple orders of magnitude.
import sys
import os
# setting path
sys.path.append('/work/scripts')
import pandas as pd
import seaborn as sns
import output_utils as ou
import xarray as xr
from glob import glob
import numpy as np
from sklearn.metrics import r2_score
from matplotlib import pyplot as plt
Specify x and y coordinate of NetCDF output
cell_y_coord = 0
cell_x_coord = 1
Soil variables are defined by layer, and soil layers are defined by the DVM-DOS-TEM layering scheme and the evolution of the soil column over the simulation. Therefore, the soil variable outputs likely will not correspond in depth with available measurements. We will need to interpolate variables output by soil layer to estimate their values at the discrete depths of available soil measurements. The following function is based on a script by Helene Genet.
depthlist = [0.05, 0.1, 0.2, 0.3] #define depths to interpolate, based on measurement depths available
def get_lwclayer_tlayer(depthlist, run_dir, var):
#interpolate output by layer to specified list of depths
#args:
# depthlist (List[float]): list of depths to interpolate variables to, must be positive and > 0
# run_dir (str): path to parent directory of run
# var (str): name of output variable output by layer (e.g. 'LWCLAYER', 'TLAYER'). Function currently supports monthly transient output
#return pd.DataFrame() with columns 'time', 'x', 'y', 'layer', 'z', 'type', '[var]'
### read the netcdf output files and compute year from the time dimension
data = xr.open_dataset(f'{run_dir}/output/{var}_monthly_tr.nc')
data = data.to_dataframe()
data.reset_index(inplace=True)
data.dtypes
data['time'] = data['time'].astype('|S80')
data['time'] = data['time'].astype('|datetime64[ns]')
data['month'] = data['time'].dt.month
data['year'] = data['time'].dt.year
data = data.sort_values(['time','x','y','layer'])
### read the netcdf output files on soil structure and compute year from the time dimension
dz = xr.open_dataset(f'{run_dir}/output/LAYERDZ_monthly_tr.nc')
dz = dz.to_dataframe()
dz.reset_index(inplace=True)
dz.dtypes
dz['time'] = dz['time'].astype('|S80')
dz['time'] = dz['time'].astype('|datetime64[ns]')
dz['month'] = dz['time'].dt.month
dz['year'] = dz['time'].dt.year
dz = dz.sort_values(['time','x','y','layer'])
### read the netcdf output files on soil structure and compute year from the time dimension
lt = xr.open_dataset(f'{run_dir}/output/LAYERTYPE_monthly_tr.nc')
lt = lt.to_dataframe()
lt.reset_index(inplace=True)
lt.dtypes
lt['time'] = lt['time'].astype('|S80')
lt['time'] = lt['time'].astype('|datetime64[ns]')
lt['month'] = lt['time'].dt.month
lt['year'] = lt['time'].dt.year
lt = lt.sort_values(['time','x','y','layer'])
dz=pd.merge(dz, lt[['LAYERTYPE', 'time', 'x', 'y', 'layer']], on=['time','x','y','layer'])
### compute the depth of the bottom of every layers
dz['z'] = dz.groupby(['time','x','y'])['LAYERDZ'].cumsum(axis=0)
### loop through the list of depths of reference to compute the soil variable at that depth via linear interpolation
stdz = []
for i in range(len(depthlist)):
dpth = depthlist[i]
print("depth:", dpth,"m")
# extract the top and bottom layers the closest to the depth of reference
dz['diff'] = dz['z']-float(dpth)
top = dz.loc[dz[(dz['diff'] <= 0)].groupby(['time','x','y'])['diff'].idxmax()]
bot = dz.loc[dz[(dz['diff'] >= 0)].groupby(['time','x','y'])['diff'].idxmin()]
# select the variable value for each of these top and bottom layers
datatop = pd.merge(data, top[['year','month', 'x','y','layer','LAYERDZ','LAYERTYPE','z']], how="left", on=['layer','year','month', 'x','y'])
datatop = datatop[datatop['z'].notna()]
datatop = datatop.rename(columns={"layer": "layertop", var: var+"top", "LAYERDZ": "dztop", "z": "ztop", "LAYERTYPE": "typetop"})
databot = pd.merge(data, bot[['year', 'month', 'x','y','layer','LAYERDZ','LAYERTYPE','z']], how="left", on=['layer','year','month' , 'x','y'])
databot = databot[databot['z'].notna()]
databot = databot.rename(columns={"layer": "layerbot", var: var+"bot", "LAYERDZ": "dzbot", "z": "zbot", "LAYERTYPE": "typebot"})
# merge the data to do the linear interpolation
datastdz = pd.merge(datatop, databot, how="outer", on=['time','year', 'month', 'x','y'])
datastdz['a'] = (datastdz[var+"top"] - datastdz[var+"bot"]) / (datastdz['ztop'] - datastdz['zbot'])
datastdz['b'] = datastdz[var+"top"] - (datastdz['a'] * datastdz['ztop'])
datastdz[var] = (datastdz['a'] * float(dpth)) + datastdz['b']
datastdz['z'] = float(dpth)
datastdz['layer'] = i
datastdz['type'] = datastdz['typebot']
datastdz = datastdz[['time','x','y','layer','z','type',var]]
stdz.append(datastdz)
stdz = pd.concat(stdz)
return stdz
Create list of working directories based on sensitivity analysis parent directory.
out_dir ='/data/workflows/US-Prr_SWC_SA/'
run_dirs = [d for d in glob(out_dir+'*/', recursive = True) if 'sample' in d]
Loop through SA runs and create dataframes for LWCLAYER, TLAYER, and other outputs of interest (GPP, RH, etc.)
As we can see, carbon fluxes are highly sensitive to the soil thermal/hydrological parameters.
Parameter Sensitivity
To better understand the sensitivity of DVM-DOS-TEM output variables to changing the thermal/hydrological parameters, we can plot the relationships in a grid.
First load in the sample matrix.
reacted with thumbs up emoji reacted with thumbs down emoji reacted with laugh emoji reacted with hooray emoji reacted with confused emoji reacted with heart emoji reacted with rocket emoji reacted with eyes emoji
-
Calibration of a new community type within DVM-DOS-TEM often involves tuning vegetation and soil carbon related parameters to match observations or literature compilations of carbon stocks and fluxes. While such calibration is necessary to begin representing ecosystem dynamics, it is still only part of the picture. In DVM-DOS-TEM, the thermal and hydrological properties of the soil exert considerable influence on carbon stocks and fluxes through processes affected by soil moisture and temperature. It is recommended that the thermo-hydrological parameter space is explored during calibration, to best represent important processes affecting ecosystem dynamics.
This example provides a framework for running sensitivity analysis on parameters with the greatest influence on the thermal and hydrological behavior of the model. In this case, we will be comparing model output soil moisture and temperature with site-based observations, in an attempt to determine an optimal set of parameters.
A list of important parameters with units, ranges, and literature references was developed by Helene Genet.
Parameters influencing soil hydro/thermal state
cmt_envground.txt
tcsolid(m): Soil thermal conductivity for moss (W/mK)
current value: 0.25, range: 0.005-0.5
Ekici et al. 2015,Jiang et al. 2015, O’Donnell et al. 2009
tcsolid(f): Soil thermal conductivity for fibric (W/mK)
current value: 0.25, range: 0.005-0.5
Ekici et al. 2015,Jiang et al. 2015, O’Donnell et al. 2009
tcsolid(h): Soil thermal conductivity for humic (W/mK)
current value: 0.25, range: 0.02-2.0
Ekici et al. 2015,Jiang et al. 2015, O’Donnell et al. 2009
porosity(m): porosity for moss layers (m3/m3)
current value: 0.98, range: 0.85-0.99
O’Donnell et al. 2009
porosity(f): porosity for fibric layers (m3/m3)
current value: 0.95, range: 0.85-0.99
O’Donnell et al. 2009
porosity(h): porosity for humic layers (m3/m3)
current value: 0.8, range: 0.7-0.9
O’Donnell et al. 2009
bulkden(m): bulk density for moss (g/m3)
current value: 25,000, range: 10,000 - 80,000
Tuomi et al. 2020, Rodionov et al. 2007, O’Donnell et al. 2009
bulkden(f): bulk density for fibric (g/m3)
current value: 51,000, range: 20,000 - 200,000
Tuomi et al. 2020, Rodionov et al. 2007, O’Donnell et al. 2009
bulkden(h): bulk density for humic (g/m3)
current value: 176,000, range: 100,000 - 800,000
Tuomi et al. 2020, Rodionov et al. 2007, O’Donnell et al. 2009
hksat(m): hydraulic conductivity at saturation for moss (mm/s)
current value: 0.15, range: 0.0002 - 30
Ekici et al. 2015, Letts et al. 2000, Liu et al. 2019
hksat(f): hydraulic conductivity at saturation for fibric (mm/s)
current value: 0.28, range: 0.0002 - 30
Ekici et al. 2015, Letts et al. 2000, Liu et al. 2019
hksat(h): hydraulic conductivity at saturation for humic (mm/s)
current value: 0.002, range: 0.00004 - 2.01
Ekici et al. 2015, Letts et al. 2000, Liu et al. 2019
nfactor(s): Summer nfactor
current value: 1.5, range: 0.2-2.0
Kade et al. 2006, Klene et al. 2001
snwalbmax
current value: 0.8, range: 0.7-0.85
Te Beest et al. 2016, Petzold et al. 1975, Loranty et al. 2011
snwalbmin
current value: 0.4, range: 0.4-0.6
Te Beest et al. 2016, Petzold et al. 1975, Loranty et al. 2011
cmt_dimground.txt
snwdenmax (kg/m3)
current value: 250, range: 100-800
Domine et al. 2016, Gerland et al. 1999, Muskett 2012
snwdennew (kg/m3)
current value: 50, range: 10-250
Domine et al. 2016, Gerland et al. 1999, Muskett 2012
cmt_bgcsoil.txt
rhq10
current value: 2, range: 1.6-2.4
rhq10_w
current value: 2, range: 0.85-0.99
SA setup
We will use a few useful imports from the 'scripts' directory in addition to a few other packkages.
Define working directory and run setup string for the sensitivity driver.
Set up the sensitivity driver.
Parameters to adjust:
The sensitivity driver has the option to define a percent to perturb the parameters, however this does not work so well for parameters spanning multiple orders of magnitude, so we will manually specify the bounds. Additionally, we will specify which parameters should be sampled on a log-uniform scale; i.e. the parameters spanning multiple orders of magnitude.
Now we will define which outputs to activate and design the experiment.
The following code block will clean the outputs, setup the sensitivity analysis to run in parallel, and run the experiment.
Sensitivity visualization
Imports
Specify x and y coordinate of NetCDF output
Soil variables are defined by layer, and soil layers are defined by the DVM-DOS-TEM layering scheme and the evolution of the soil column over the simulation. Therefore, the soil variable outputs likely will not correspond in depth with available measurements. We will need to interpolate variables output by soil layer to estimate their values at the discrete depths of available soil measurements. The following function is based on a script by Helene Genet.
Create list of working directories based on sensitivity analysis parent directory.
Loop through SA runs and create dataframes for LWCLAYER, TLAYER, and other outputs of interest (GPP, RH, etc.)
Load in observation data and format.
Filter SA outputs to align with observations.
Merge tlayer and lwclayer dataframes to only contain time periods with data overlapping observations.
Now we can calculate r2 between observed and modeled soil temperature.
Now we can plot out the results to visualize the observations, range on simulations, and simulation with the highest r2
TLAYER
LWCLAYER
GPP
Here we are plotting GPP, highlighting the set of parameters optimized for soil temperature in black
We can also plot GPP while highlighting the set of parameters optimized for soil moisture in black.
As we can see, carbon fluxes are highly sensitive to the soil thermal/hydrological parameters.
Parameter Sensitivity
To better understand the sensitivity of DVM-DOS-TEM output variables to changing the thermal/hydrological parameters, we can plot the relationships in a grid.
First load in the sample matrix.
Next we will separate visualizations by season, so we need to extract the month from results dataframes.
Now merge the results with the sample matrix.
Summer Sensitivity
Winter Sensitivity
Beta Was this translation helpful? Give feedback.
All reactions