%% Model: Duke Forest CO2 dynamics % % This script numerically integrates the carbon dioxide balance with parameters % selected for the Milwaukee CO2 data. The following figures are produced: % % Figure 3: Comparison of measured and modeled data clear all; close all; clc %% Load data x=tdfread('DukeForest_soil_data.txt'); VWC=x.VWC; T=x.Temp_0x28C0x29; CO2=x.Cm_0x28ppm0x2F100x29/1e4*10; %convert ppm/10 to percent %% Define parameters N=length(VWC); %simulation time, 30-minute increments dt=1; %time-step, 30-minute increments n=0.55; %porosity zr=70/1000; %root depth, m L=0.05*zr; %diffusion length, m D0=0.157/100^2/L*3600/2; %Diffusion coefficient, m^2/30-min r=0.88; %respiration rate, umol/mol C/30-min Q10=exp(1); %Q10 coefficient, [] ks=-0.5; %soil moisture coefficient, [] CO2atm=410/1e4; %Atmospheric concentration (%) kH=0.8317; %Henry's Law Constant, [] R=8.3145; %Ideal gas constant, J/kg/mol s=VWC/n; %convert VWC to s %% Calculate adjustment functions, retardation factor, and timescales %Diffusion fDs=n^2.5*(1-s).^1.5; fDT=((T+273)/293).^1.75; %Respiration fRT=Q10.^((T-25)/10); fRs=exp(ks*s); %Retardation factor (assumed constant equal to the mean) Rf=((1-s)+kH*s*R.*T)*n*zr; Rf=mean(Rf)*ones(1,length(s)); %Linear model timescales taud=D0/n/zr/mean(Rf); taur=r/n/zr/mean(Rf); %% Model integration c=zeros(1,N-1); c(1)=CO2(1); CO2bar=mean(CO2); for i=1:N-1 %Diffusion flux Fdiff(i)=fDs(i)*fDT(i)*D0*(c(i)-CO2atm); %Respiration flux Fresp(i)=fRs(i)*fRT(i)*r; %Forward Euler integration c(i+1)=c(i)+dt*(-Fdiff(i)+Fresp(i))/Rf(i); end save co2model.mat %% Plot time-series figures for paper figure(5) set(gcf,'Units','inches','PaperPosition',[0 0 3 3]) set(gcf,'Color',[1 1 1]) plot((1/48):(1/48):(length(CO2)/48),CO2,'o','Color',[0.8 0.8 0.8]) hold all plot((1/48):(1/48):(length(CO2)/48),c,'k-') set(gca,'FontSize',20) xlabel('Time (d)') ylabel('CO_2 (%)') xlim([400 450]) xticklabels(0:10:50) ylim([0 1.5]) legend('Measured','Modeled','location','southwest') title('(c)') figure(6) set(gcf,'Units','inches','PaperPosition',[0 0 3 3]) set(gcf,'Color',[1 1 1]) plot(CO2,c,'o','Color',[0.8 0.8 0.8]) hold all plot([0 1.6],[0 1.6],'k--') set(gca,'FontSize',20) xlim([0 1.6]) yticks(0:0.5:1.5) ylim([0 1.6]) xlabel('Measured CO_2 (%)') ylabel('Modeled CO_2 (%)') title('(d)')