options(warn=-1)
# Import dependencies
library(rgdal)
library(raster)
library(gdalUtils)
library("readxl")
library(viridis)
library(plyr)
library(sp)
library(ggplot2)
library(rasterVis)
library(gridExtra)
library(grid)
# Define function to adjust two rasters
adj_raster = function(x,y,crop_extent=NULL){
# Decide which raster has lower resolution
if(res(x)[1] < res(y)[1]){
low_res = y
high_res = x
} else {
low_res = x
high_res = y
}
if(is.null(crop_extent)){
# Calculate the shared extent between the rasters
crop_ext = intersect(extent(low_res),extent(high_res))
# Crop both rasters to the shared extent
crop_low_res = crop(low_res,crop_ext)
crop_high_res = crop(high_res,crop_ext)
} else {
crop_low_res = crop(low_res,crop_extent)
crop_high_res = crop(high_res,crop_extent)
}
# Create a new raster to resample the high resolution raster into
resample_high_res = raster(nrow=nrow(crop_high_res)*(res(crop_high_res)/res(crop_low_res))[1],
ncol=ncol(crop_high_res)*(res(crop_high_res)/res(crop_low_res))[1],
xmn=extent(crop_high_res)[1],
xmx=extent(crop_high_res)[2],
ymn=extent(crop_high_res)[3],
ymx=extent(crop_high_res)[4])
# Resample the high resolution raster
resample_high_res = resample(crop_high_res,resample_high_res)
# Return results in the correct order
if(res(x)[1] < res(y)[1]){
return(list(resample_high_res,crop_low_res))
} else {
return(list(crop_low_res,resample_high_res))
}
}
# Define a function that calculates the geometric mean
gmean = function(x,y){
return(exp(mean(log(c(x,y)))))
}
# Define a function to calculate the multiplicative 95% confidence interval of a list of measurements
mul_CI = function(vals){
sem = 10^(sd(log10(vals))/sqrt(length(vals))*1.96)
std = 10^(sd(log10(vals))*1.96)
return(gmean(sem,std))
}
# Define a function to propagate uncertainties through a product
CI_prod_prop = function(mul_CIs){
return(10^sqrt(sum(log10(mul_CIs)^2)))
}
# Define a funtion to overlay data in maps that may be in differnt sizes
overlay_maps = function(x,y, calc_sum=FALSE,crop_extent=NULL){
# Adjust the two rasters
tmp = adj_raster(x,y,crop_extent=crop_extent)
res_x = tmp[[1]]
res_y = tmp[[2]]
# Multiply the values in the two maps
overlayed_map = overlay(res_x,res_y,fun=function(x,y){x*y})
if(calc_sum){
moll = '+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs'
# Project to an equal area projection
overlayed_moll = projectRaster(overlayed_map,crs=moll)
#Calculate the total leaf mass
tot_sum = cellStats(overlayed_moll, stat = 'sum')*prod(res(overlayed_moll))
return(tot_sum)
}else {
return(overlayed_map)
}
}
In order to estimate the total leaf area we rely on two procedures. The first is based on remotely sensed measurement of the leaf area index (LAI), which measures the total area of leaves per unit surface area. The second is based on ground based measurements of LAI.
We rely on the GLASS Leaf Area Index (LAI) product for this estimate (http://glcf.umd.edu/data/lai/ ). We use 8-day composites taken throughout 2015 and calculate the total leaf area per sample:
# Create results file
result_MODIS = data.frame(file_name=character(), day=numeric(), lai = numeric())
# CRS with equal area
moll = '+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs'
# Try loading the results of the analysis. If The results are not found, run the analysis
tryCatch({result_MODIS = read.csv('../data/MODIS_LAI_sum.csv'); }, error = function(err){
# Find all HDF files in the MODIS dataset
files = list.files('../data/GLASS_LAI_data/MODIS/')
for(file in files){
print(file)
# Load raster
glass_data = get_subdatasets(paste('data/GLASS_LAI_data/MODIS/',file,sep=''))
lai = readGDAL(glass_data[1])
lai_raster = raster(lai)
# Reproject the raster to equal area projection
molllai <- projectRaster(lai_raster, crs=moll)
# Calculate the total leaf area
total_leaf_area = cellStats(molllai,stat = 'sum')*prod(res(molllai))
# Store results
result_MODIS = rbind(result_MODIS,data.frame(file_name=file,day=substr(unlist(strsplit(file,'\\.'))[3],6,10),lai=total_leaf_area))
}
# Save results
write.csv(result_MODIS,'../data/MODIS_LAI_sum.csv')
}
)
# Find the raster with a total area closest to the annual mean
best_file = result_MODIS[which.min(abs(result_MODIS$lai-mean(result_MODIS$lai))),2]
sprintf('We estimate the total leaf area based on GLASS LAI product at: %.1e m^2',mean(result_MODIS$lai))
# Open the raster with total leaf area closest to the mean
best_lai_raster = raster(get_subdatasets(paste('../data/GLASS_LAI_data/MODIS/',best_file,sep=''))[1])
Our second method for estimating the total leaf area is based on field measurement of LAI in different biomes. This data takes into account only places which have vegetation in them, and thus we need aditional data on the percent coverage of surface area by vegetation in each biome. For the field measured LAI data we use Asner et al..
asner_LAI_data = read_excel('../data/literature_data.xlsx','MODIS_land_use',skip = 1)
asner_LAI_data
We use biome data from MODIS land cover product (http://glcf.umd.edu/data/lc/ ). For each biome, we replace all the pixels whithin the biome with the mean LAI measured in the biome based on the table produced by Asner et al.
# Load biome map raster
biome_data = raster('../data/LC_5min_global_2012.tif')
# Map the biome codes in the map to the mean LAI measured by Asner et al.
asner_lai_map = mapvalues(as.matrix(biome_data),from = asner_LAI_data$Value, to = asner_LAI_data$`LAI [m^2 m^-2]`)
Not all locations on in each biome contain vegetation. To take this factor into account, we use data on the tree cover and short vegetation conver in each location in the world, taken from Song et al..
# Load tree cover data
TG_data = raster('../data/VCF5KYR_2016001_001_2018224210310.tif',band=1)/100
# Load short vegetation cover data
SV_data = raster('../data/VCF5KYR_2016001_001_2018224210310.tif',band=2)/100
# Calculate the fraction of land covered by vegetation
plant_frac_data = (SV_data+TG_data)
asner_lai_raster = raster(asner_lai_map)
# Convert data to raster
extent(asner_lai_raster) = extent(biome_data)
res(asner_lai_raster) = res(biome_data)
crs(asner_lai_raster) = crs(biome_data)
# Multiply the map of the fraction of vegetation by the map of the ground based LAI measurements
asner_leaf_area_raster = overlay_maps(plant_frac_data,asner_lai_raster)
# Project to an equal area projection
asner_leaf_area_raster_moll = projectRaster(asner_leaf_area_raster,crs=moll)
#Calculate the total leaf area
asner_total_leaf_area = cellStats(asner_leaf_area_raster_moll, stat = 'sum')*prod(res(asner_leaf_area_raster_moll))
sprintf('We estimate the total leaf area based on field measurements at: %.1e m^2',asner_total_leaf_area)
As our best estimate for the total leaf area, we use the geometric mean of the estimates based on remote sensing and field measurements:
# Calculate the geometric mean of the estimates based on remote sensing nad ground measurements
best_global_leaf_area = gmean(asner_total_leaf_area,mean(result_MODIS$lai))
sprintf('Our best estimate for the total leaf area is %.1e m^2',best_global_leaf_area)
In order to convert leaf area to leaf mass, we rely on two methods. The first is based on the Glopnet database (http://bio.mq.edu.au/~iwright/glopian.htm). We use the geometric mean of leaf mass per area (LMA) in g per m^2, and multiply it by the our best estimate for the total leaf area
# Load glopnet data
glopnet_data = read_excel('../data/literature_data.xlsx','glopnet_data',skip = 1)
# calculate the geometric mean of the LMA
glopnet_mean = 10^mean((glopnet_data$`log LMA`), na.rm = TRUE)
# Multiply our best estimate of the total leaf area by the mean LMA to get the total mass of leaves
method1_mean = best_global_leaf_area*glopnet_mean
sprintf('Our estimate for the total leaf mass based on glopnet is %.0f Gt',method1_mean/1e15)
The second method we use is a recent global mapping of plant traits Butler et al.. We use a map of the specific leaf area (SLA), which is measured in m$^2$ per kg of dry leaf mass. We convert the SLA into units of gram dry leaf mass per m$^2$:
# Load SMA data
SLA_data = read.csv('../data/spat_1_superpft_sla_large.csv')
# Convert units from m^2 per kg to g per m^2
SLA_data$mean = 1000/SLA_data$mean
# Convert the coordinate data into coordinates
SLA_df = data.frame(lat=SLA_data$lat, lon=SLA_data$lon,mean=SLA_data$mean)
coordinates(SLA_df) <- ~lon+lat
# Set the projection of the SLA map to lat-lon
proj4string(SLA_df)=CRS("+init=epsg:4326")
# Create a raster of SLA
gridded(SLA_df) = TRUE
SLA_raster = raster(SLA_df)
We apply the map of leaf mass per area on both the maps of leaf area we constructed in the previous section (based on remote sensing or ground measuremetns)
# Calculate the leaf mass based on ground measurements of LAI
# Multiply the mass of the leaf mass per unit area by the map of leaf area
asner_leaf_mass_raster = overlay_maps(asner_leaf_area_raster,SLA_raster)
# Project to an equal area projection
asner_leaf_mass_raster_moll = projectRaster(asner_leaf_mass_raster,crs=moll)
#Calculate the total leaf area
asner_total_leaf_mass = cellStats(asner_leaf_mass_raster_moll, stat = 'sum')*prod(res(asner_leaf_mass_raster_moll))
# Calculate the leaf mass based on remote sensing measurements of LAI
# Multiply the mass of the leaf mass per unit area by the map of leaf area
glass_leaf_mass_raster = overlay_maps(best_lai_raster,SLA_raster)
# Project to an equal area projection
glass_leaf_mass_raster_moll = projectRaster(glass_leaf_mass_raster,crs=moll)
#Calculate the total leaf mass
glass_total_leaf_mass = cellStats(glass_leaf_mass_raster_moll, stat = 'sum')*prod(res(glass_leaf_mass_raster_moll))
# Calculate the geometric mean of the two estimates for the global mass of leaves
method2_mean = gmean(glass_total_leaf_mass,asner_total_leaf_mass)
sprintf('Our best estimate for the total mass of leaves based on a map of leaf mass per unit area is %.0f Gt', method2_mean/1e15)
As our best estimate for the global leaf mass based on remote sensing, we use the geometric mean of the our estimates based on values from the glopenet database and the map of the specific leaf area
# Calculate the geometric mean of the estimates for leaf mass based on glopnet and the map of leaf mass
best_leaf_mass_remote_sensing = gmean(method1_mean,method2_mean)
sprintf('Our best estimate for the total mass of leaves is %.0f Gt', best_leaf_mass_remote_sensing/1e15)
In this section, we estimate the fraction of the global leaf mass that is found in woody plants and herbaceous plants. In order to generate the estimate, we use four different methodologies - with two different leaf area maps (the ground-based LAI map and the remote sensing-based LAI map) and two different conversions from leaf area to leaf mass (one based on glopnet and one based on plant trait maps in Butler et al.).
# Produce a mask with the surface area that is covered by herbaceous plants
# (grasslands, croplands and cropland/Natural vegetation mosaics)
herb_raster = (biome_data == 10) | (biome_data == 12) | (biome_data == 14)
## Method 1
# Estimate the fraction of herbaceous leaf mass using ground-based leaf area map and
# a constant leaf mass per unit leaf area based on glopnet
#Overlay the maps of herbaceous cover with the ground-based leaf area map
ground_constant_herb_mass = overlay_maps(herb_raster,asner_lai_raster, calc_sum=TRUE)
# Calculate the fraction of leaf mass in herbaceous plants
ground_constant_herb_mass_frac = ground_constant_herb_mass/asner_total_leaf_area
## Method 2
# Estimate the fraction of herbaceous leaf mass using remote sensing-based leaf area map and
# a constant leaf mass per unit leaf area based on glopnet
#Overlay the maps of herbaceous cover with the remote sensing-based leaf area map
rs_constant_herb_mass = overlay_maps(herb_raster,best_lai_raster, calc_sum=TRUE)
# Calculate the fraction of leaf mass in herbaceous plants
rs_constant_herb_mass_frac = rs_constant_herb_mass/mean(result_MODIS$lai)
## Method 3
# Estimate the fraction of herbaceous leaf mass using ground-based leaf mass map based
# on a the leaf mass per unit area map
#Overlay the maps of herbaceous cover with the ground-based leaf mass map
ground_butler_herb_mass = overlay_maps(herb_raster,asner_leaf_mass_raster, calc_sum=TRUE)
# Calculate the fraction of leaf mass in herbaceous plants
ground_butler_herb_mass_frac = ground_butler_herb_mass/asner_total_leaf_mass
## Method 4
# Estimate the fraction of herbaceous leaf mass using remote sensing-based leaf mass
# map based on a the leaf mass per unit area map
#Overlay the maps of herbaceous cover with the remote sensing-based leaf mass map
rs_butler_herb_mass = overlay_maps(herb_raster,glass_leaf_mass_raster, calc_sum=TRUE)
# Calculate the fraction of leaf mass in herbaceous plants
rs_butler_herb_mass_frac = rs_butler_herb_mass/glass_total_leaf_mass
# Calculate the geometric mean of all methods
best_herb_mass_frac = prod(c(ground_constant_herb_mass_frac,rs_constant_herb_mass_frac,ground_butler_herb_mass_frac,rs_butler_herb_mass_frac))^(1/4)
sprintf('Our best estimate fraction of leaf mass that is in herbaceous plants is %.0f percent',best_herb_mass_frac*100)
To estimate the total mass of leaves of C4 plants, we estimate seperately the mass leaves of C4 crops and C4 natural vegetation
To estimate the total mass of leaves of C4 crops, we rely on the leaf masss maps we constructed in previous sections, in conjunction with maps of the fraction of area that contains C4 crops from Monfreda et al.. We use the dominant C4 crops maize, sorghum and sugarcane as an indicator for the amount of C4 crops at every location.
# Load data on C4 crops
maize = raster('../data/C4_leaf_mass/HarvestedAreaYieldMajorCrops_Geotiff/maize_HarvAreaYield_Geotiff/maize_HarvestedAreaFraction.tif')
sorghum = raster('../data/C4_leaf_mass/HarvestedAreaYieldMajorCrops_Geotiff/sorghum_HarvAreaYield_Geotiff/sorghum_HarvestedAreaFraction.tif')
sugarcane = raster('../data/C4_leaf_mass/HarvestedAreaYieldMajorCrops_Geotiff/sugarcane_HarvAreaYield_Geotiff/sugarcane_HarvestedAreaFraction.tif')
C4_crops = (maize+sorghum+sugarcane)
To generate our estimate of the total mass of leaves of C4 crops, we multiply the map of the fraction of surface area covered by C4 crops by our four maps of the total mass of leaves
total_C4_crop_leaf_mass1 = overlay_maps(C4_crops,asner_lai_raster, calc_sum = TRUE)*100
total_C4_crop_leaf_mass2 = overlay_maps(C4_crops,best_lai_raster, calc_sum = TRUE)*100
total_C4_crop_leaf_mass3 = overlay_maps(C4_crops,glass_leaf_mass_raster, calc_sum = TRUE)
total_C4_crop_leaf_mass4 = overlay_maps(C4_crops,asner_leaf_mass_raster, calc_sum = TRUE)
#glass_total_C4_leaf_mass/rs_butler_herb_mass
best_tot_C4_crop_leaf_mass = prod(c(total_C4_crop_leaf_mass1,total_C4_crop_leaf_mass2,total_C4_crop_leaf_mass3,total_C4_crop_leaf_mass4))^(1/4)
Our best estimate of the total mass of leaves in C4 crops is the geometric mean of the four different estimates:
best_tot_C4_crop_leaf_mass = prod(c(total_C4_crop_leaf_mass1,total_C4_crop_leaf_mass2,total_C4_crop_leaf_mass3,total_C4_crop_leaf_mass4))^(1/4)
sprintf('Our best esimate for the total mass of leaves of C4 crops is ≈%.1f Gt', best_tot_C4_crop_leaf_mass/1e15 )
To estimate the total mass of leaves of natural C4 plants, we rely on a depiction of the fraction of C4 coverage globally (Still et al.). Still et al. includes C4 coverage of biomes that are dominated by woody platns such as savanna and shrublands, because these environments include also herbaceous understory. Currently, there is no procedure to directly estimate the fraction of the understory out of the leaf mass in these biomes. Therefore, we generate two estimates of the fraction of C4 plants out of the global leaf mass - one that takes into account the fraction of C4 plants only out of grasslands, and the other includes also savanna and shrublands.
We generate the a map of C4 area coverage only in grasslands. We use the natural C4 plant coverage map in conjunction with the leaf masss maps we constructed in previous sections to generate an estiamte for the total leaf mass of natural C4 plants.
# Load the data from Still et al.
natural_c4_frac = read.asciigrid('../data/C4_leaf_mass/ISLSCP_C4_1DEG_932/data/c4_percent_1d.asc')
natural_c4_frac = raster(natural_c4_frac)/100
natural_c4_frac[natural_c4_frac<0] = 0
crs(natural_c4_frac) = CRS("+proj=longlat +datum=WGS84")
grasslands_raster = biome_data
grasslands_raster = grasslands_raster==10
corrected_natural_c4_frac_grassland = overlay_maps(natural_c4_frac,grasslands_raster)
nat_c4_mass_frac_grassland1 = overlay_maps(corrected_natural_c4_frac_grassland,asner_lai_raster,calc_sum = TRUE)*100
nat_c4_mass_frac_grassland2 = overlay_maps(corrected_natural_c4_frac_grassland,best_lai_raster,calc_sum = TRUE)*100
nat_c4_mass_frac_grassland3 = overlay_maps(corrected_natural_c4_frac_grassland,asner_leaf_mass_raster,calc_sum = TRUE,crop_extent = extent(c(-180,180,-56,83)))
nat_c4_mass_frac_grassland4 = overlay_maps(corrected_natural_c4_frac_grassland,glass_leaf_mass_raster,calc_sum = TRUE,crop_extent = extent(c(-180,180,-56,83)))
best_nat_c4_mass_frac_grassland = prod(nat_c4_mass_frac_grassland1,nat_c4_mass_frac_grassland2,nat_c4_mass_frac_grassland3,nat_c4_mass_frac_grassland4)^(1/4)
We generate the a map of C4 area coverage in grasslands and also savanna and shrublands. We use the natural C4 plant coverage map in conjunction with the leaf masss maps we constructed in previous sections to generate an estiamte for the total leaf mass of natural C4 plants.
# Create a mask of land excluding croplands
crop_raster = biome_data
crop_raster[crop_raster==0] = 12
crop_raster[crop_raster!=12] = 1
crop_raster[crop_raster==12] = 0
corrected_natural_c4_frac = overlay_maps(natural_c4_frac,crop_raster)
nat_c4_mass_frac_all1 = overlay_maps(corrected_natural_c4_frac,asner_lai_raster,calc_sum = TRUE)*100
nat_c4_mass_frac_all2 = overlay_maps(corrected_natural_c4_frac,best_lai_raster,calc_sum = TRUE)*100
nat_c4_mass_frac_all3 = overlay_maps(corrected_natural_c4_frac,asner_leaf_mass_raster,calc_sum = TRUE,crop_extent = extent(c(-180,180,-56,83)))
nat_c4_mass_frac_all4 = overlay_maps(corrected_natural_c4_frac,glass_leaf_mass_raster,calc_sum = TRUE,crop_extent = extent(c(-180,180,-56,83)))
best_nat_c4_mass_frac_all = prod(nat_c4_mass_frac_all1,nat_c4_mass_frac_all2,nat_c4_mass_frac_all3,nat_c4_mass_frac_all4)^(1/4)
Our best estimate for the mass of leaves in C4 plants in natural habitats is the geometric mean of the two estimates using only grasslands or all natural biomes:
best_nat_c4_leaf_mass = gmean(best_nat_c4_mass_frac_all,best_nat_c4_mass_frac_grassland)
sprintf('Our best esimate for the total mass of leaves of natural C4 plants is ≈%.1f Gt', best_nat_c4_leaf_mass/1e15 )
As our best estimate of the total mass of leaves of C4 plants, we sum the our best estimates for the total mass of leaves of C4 crops and natural C4 plants
best_tot_C4_leaf_mass = best_nat_c4_leaf_mass+best_tot_C4_crop_leaf_mass
sprintf('Our best esimate for the total mass of leaves of C4 plants is ≈%.1f Gt', best_tot_C4_leaf_mass/1e15 )
We also estimate the fraction of the total herbaceous leaf mass that is contributed by C4 plants:
tot_herb_leaf_mass = prod(ground_constant_herb_mass*100,rs_constant_herb_mass*100,ground_butler_herb_mass,rs_butler_herb_mass)^(1/4)
best_tot_C4_leaf_frac = best_tot_C4_leaf_mass/tot_herb_leaf_mass
sprintf('Our best esimate for the fraction of total mass of leaves of herbaceous plants that is in C4 plants is ≈%.0f percent', best_tot_C4_leaf_frac*100 )
We calculate an uncetainty projection for our estimate of the total mass of leaves based on remote sensing. This is a multiplicative uncertainty akin to a 95% confidence interval. We start with estimating the uncertainty around our estimate of the total area of leaves:
leaf_area_CI = mul_CI(c(asner_total_leaf_area,mean(result_MODIS$lai)))
sprintf('Our best projection for the uncertainty associated with our estimate of the total leaf area is %.1f-fold', leaf_area_CI)
In order to calculate the total mass of leaves based on remote sensing, we multiply the total leaf mass by the leaf mass per unit leaf area. Next, we calculate the uncertainty associated with our estimate of the leaf mass per unit leaf area. To calculate the characteristic leaf mass per unit leaf area, we relied on an estimate based on a database of plant traits in ≈2000 plant species, and a map of plant traits. For each estimate, we calculate the uncertainty associated with it.
Our estimate was based on the geomteric mean across plant species. As an estimate of the uncertainty associated with this estimate, we calculate the 95% confidence interval of around the geometric mean of the leaf mass per unit leaf area across all species:
glopnet_CI = mul_CI(na.exclude(glopnet_data$`log LMA`))
sprintf('Our best projection for the uncertainty associated with our estimate of the leaf mass per unit leaf area based on the species database is %.1f-fold', glopnet_CI)
The uncertainty of 1.2-fold in the estimate of the leaf mass is larger than the uncertainty between our different methodologies (the species based database and the trait maps). Therefore, we use this uncertainty as our best projection of the uncertainty associated with the characteristic mass of leaves per unit leaf area.
Combining the uncertainty associated with the total area of leaves with the uncertainty associated with our estimate for leaf mass per unit leaf area, we get a projection for the uncertainty associated with the estimate of the total mass of leaves based on remote sensing:
rs_leaf_mass_CI = CI_prod_prop(c(glopnet_CI,leaf_area_CI))
sprintf('Our best projection for the uncertainty associated with our estimate of the total leaf mass based on remote sensing is %.1f-fold', rs_leaf_mass_CI)
In order to evaluate the uncertainty associated with our estimate of the total mass of leaves of herbaceous plants, we calculate the 95% confidence interval of around the geometric mean of the mass of herbaceous leaves out of our four estimates:
herb_leaf_mass_CI = mul_CI(c(ground_constant_herb_mass*100,rs_constant_herb_mass*100,ground_butler_herb_mass,rs_butler_herb_mass))
sprintf('Our best projection for the uncertainty associated with our estimate of the total herbaceous leaf mass is %.1f-fold', herb_leaf_mass_CI)
In order to evaluate the uncertainty associated with our estimate of the total mass of leaves of C4 plants, we calculate the 95% confidence interval of around the geometric mean of the mass of C$ leaves out of our four estimates:
C4_leaf_mass_CI = mul_CI(c(nat_c4_mass_frac_all1,nat_c4_mass_frac_all2,nat_c4_mass_frac_all3,nat_c4_mass_frac_all4))
sprintf('Our best projection for the uncertainty associated with our estimate of the total C4 leaf mass is %.1f-fold', C4_leaf_mass_CI)
We also estimate the uncertainty associated with our estimate of the fraction of hebaceous leaves that is in C4 plants:
C4_leaf_mass_frac_CI = CI_prod_prop(c(herb_leaf_mass_CI,C4_leaf_mass_CI))
sprintf('Our best projection for the uncertainty associated with our estimate of the fraction of the total herbaceous leaf mass that is in C4 plants is %.1f-fold', C4_leaf_mass_frac_CI)