co2_soil.m

clc
close all
clearvars


%% %% Water Dynamics

%% Surface water runoff

% USDA Soil Conservation Service runoff curve number
CN = 1;

% Rainfall and/or irrigation (in meters per hour per square meter)
RI = 1;

% Eq. 2 - Soil water retention parameter
S = 1000./CN - 10;

% Eq. 1 - Surface water runoff (in cubic meters per hour per square meter)
Wrunoff = (RI - 0.2*S).^2./(RI + 0.8*S);

%% Soil water vertical infiltration

% Drainage coefficient (per hour)
alpha = 1;

% Volumetric water content (in cubic meters per cubic meter)
theta = 1;

% Field water capacity (in cubic meters per cubic meter)
fc = 1;

% Eq. 3 - Soil water vertical infiltration (in cubic meters per hour per square meter)
Lwater = alpha*(theta - fc);

%% Soil water evapotranspiration

% Soil surface evaporation rate (in cubic meters per hour per square meter)
Eevap = 1;

% Leaf water transpiration rate (in cubic meters per hour per square meter)
Etranp = 1;

% Diurnal factor
fd = 1;

% Eq. 4 - Evapotranspiration rate (in cubic meters per hour per square meter)
ET = (Eevap + Etranp).*fd;

%% %% Temperature Dynamics

% Time (in hours)
% t = linspace(0,168,10000);
t = linspace(0,10*24,1000);
% t = 0;

% Soil depth (in meters)
x = 0.3;

% Table 1 - Average soil temperature (in degrees Celsius)
Ta = -8e-7*t.^2 + 0.0079*t + 2.422;

% Amplitude of surface temperature fluctuation (in degrees Celsius)
Ao = 10 + 0*randn(size(t));

% Radial frequency
omega = 2*pi/24; % 1 period a day

% Characteristic depth or damping depth (in meters)
d = 0.15;

% Eq. 5 - Diurnal variation of soil temperature (in degrees Celsius)
Tfunc = @(x,t,Ta,Ao,omega,d) Ta + Ao.*exp(-x./d).*sin(omega.*t - x./d);

% Soil temperature for a given (x,t) (in degrees Celsius)
T = Tfunc(x,t,Ta,Ao,omega,d);

% plot(t,T)
% axis([0, 168, 0, 6])

%% %% CO2 Dynamics

%% Soil CO2 emission into the near surface atmosphere

% Reference CO2 diffusion coefficient at a reference temperature
Dr = 0.004;

% Reference soil temperature
Tr = 20;

% Eq. 7 - Diffusion coefficient of CO2 (in square meters per hour)
DCO2 = Dr.*(T./Tr).^2;

% Table 1 - Initials soil CO2 mass (in miligrams)
MCO2 = 99.3708;

% Land surface area (in square meters)
Aland = 10e3;

% Soil volume (in cubic meters)
Vsoil = Aland.*x;

% Concentration of CO2 (in milligrams per cubic meter)
CCO2soil = MCO2./Vsoil;

% Table 1 - Concentration of CO2 (in milligrams per cubic meter)
CCO2atm = 0.6134e-3;

% Soil depth interval (in centimeters) of interest.
deltaX = 100;

% Eq. 6 - Soil CO2 flux (in milligrams per hour)
JCO2 = DCO2 .* (CCO2soil - CCO2atm) ./ deltaX;

%% Soil respiration

% Fine root production (in grams per square meter per year)
FRP = 220.8;

% Microbial biomass (in milligrams C per kilogram dry soil);
MB = 144.3;

% Soil organic carbon (in percent)
SOM = 2.4;

% Soil pH
pH = 5.7;

% Eq. 8 - Soil respiration rate (in grams C per square meter per year)
Rsoil = 1266 - 0.1*FRP + 12*MB - 28.2*SOM - 96*pH;

% Eq. 9 - Diurnal temperature correction factor
ft = 0.4493*exp(0.04*T);

% % Land surface area (in square meters)
% Aland = 10e3;

% Eq. 10 - Soil CO2 production rate (in grams CO2 per year)
RCO2 = 3.667*Rsoil.*Aland.*ft;

plot(t,RCO2)