PDEPE for root diffusion in soils

Hello, I am very new to MATLAB and if I'm honest I feel like I really don't know what I'm doing, I've managed to create a basic PDEPE for simulating diffusion of CO2 from roots in the soil that I can impliment different effective diffusivity coefficients in depending on soil type, with root distribution being specified by a root decay with depth function. The main issue I am having right now is that I am needing to test my model to see the different concentration profiles it produces under different conditions, taking data from papers to base it on, iterating through each parameter at a time and seeing how it effects the concentration profile. Right now I am having particular difficulty with the root emission rates function, whilst my model uses mol m⁻³ s⁻¹ most papers report emission rates as μmol m−2 s−1. Maths and coding is not my strong point and I am looking for any pointers on how I can implement my code such that it can convert these units to those used in my code. Any help/pointers appreciated.

diffusion_roots_pdepe()

function diffusion_roots_pdepe

x = linspace(0, 1, 500); %divides soil depth of 1m into 500 points
t = linspace(0, 86400, 500); %defines time points where concentration is computed
m = 0; %1D slab (cartesian coordinates)

sol = pdepe(m, @pdefun, @icfun, @bcfun, x, t);
C = sol(:, :,1); %Solve for 1 variable across all time and space points

figure
plot(x, C(end,:), 'LineWidth', 2)
xlabel('Depth (m)')
ylabel('Concentration (mol m⁻³)')
title('Final Concentration Profile')
grid on

figure
h_line = plot(x, C(1,:), 'b-', 'LineWidth', 2);%Reference line that updates
xlabel('Depth (m)')
ylabel('Concentration (mol m^{-3})')
title('Soil gas concentration profile')
grid on

ylim([min(C(:)), max(C(:))]) %Defines Y axis limits, based on min and max C 
xlim([0, 1]) %X limits simply based on depth 
h_title = title(''); %handle that updates title at each time step 

for i = 1:length(t) %Iterates over each time step of concentration profile
    set(h_line, 'YData', C(i,:)) %Updates concentration profile with Y data at each iteration
    h_title.String = sprintf('CO2 concentration — t = %.0f s', t(i)); %Updates title for number of seconds
    drawnow
end

function [c, f, s ] = pdefun(x, t, u, dudx)

    %c = storage term
    %f = flux term
    %s = source/sink   

%EFFECTIVE AIR DIFFUSIVITY

D_air = 1e-5; %Diffusion of CO2 in free air, m² s⁻¹
theta = 0.12; %Air filled porosity of soil 
Dp_Do = theta^2; %Penman relative soil diffusivity, where theta is gas filled porosity of soil
D_eff = Dp_Do * D_air; %Effective soil diffusivity scaled by Dp_Do

%MOLDRUP DIFFUSIVITY - HEAVY CLAY SOILS

epsilon = 0.12; %Air filled porosity
phi  = 0.53; %total soil porosity 
b = 11.4; %campbell soil water retention parameter

D_eff_clay_1 = D_air * (phi^2)  * (epsilon / phi)^(2 + 3/b);

%MILLINTON-QUIRK DIFFUSIVITY - HEAVY CLAY SOIL

D_eff_clay_2 = D_air * (epsilon^2) / (phi ^ (2/3));

%ROOT DISTRIBUTION

beta = 5.961; %root decay factor defining distribution with depth
R = exp(-beta * x); %Gale Grigal formula for modelling root density with depth


%root decay factor defining distribution with depth


%ROOT EMISSION 

S0 = 1e-8; % mol m⁻³ s⁻¹ — root emission rate
S = S0*R; %Scales production by root system distribution

%MICROBIAL SINK

K_micro = 0; %Microbial consumption rate, mol m⁻³ s⁻¹ (Vmax)
Km = 1e-6 ; %Saturation point, when concentration is half Vmax, mol m⁻³
%This results in the rate of increase in consumption to slow down until
%it approaches saturation and slows down 

sink = K_micro * R * u / (Km + u); %Michaelis-Menten Kinetics term 
%for non linear consumption limited by saturation, leads to saturation of
%uptake by microbes at high concentration above Km


%PDE TERMS

c = 1; %Storage coefficient 
f = D_eff * dudx; %Fickian diffusion term scaled by effective soil diffusivity
s = S - sink; %Production - Consumption 

end

%STARTING CONDITIONS
    function U0 = icfun(x) 
        U0 = 0.016; %Initial concentration, atmospheric CO2
    end

%BOUNDARY CONDITIONS 

function [pl, ql, pr, qr] = bcfun(xl, ul, xr, ur, t)
    %pl, ql, pr, qr left and right boundary coefficients respectively

    h = 0.0005; %Boundary exchange coefficient, m s⁻¹
    C_atm = 0.016; %Atmospheric concentration, mol m⁻³

    % Left boundary Robin condition
    pl = h*(ul - C_atm); %Robin coefficient, where diffusion to atmosphere
    %is limited by h, a factor at the boundary restricting diffusion 
    ql = 1;  % Flux across the boundary
    % Right boundary Neumann condition
    pr = 0; %There is no flux across the boundary, so no flux term 
    qr = 1;%Flux term active at the boundary so C is free to change
    %There is no flux across a Neumann boundary, this can lead to build up
    %of substrate

  
end

end

%structure of code was based on matlab documentation on how to use PDEPE

Answer 1

[Torsten]

In my opinion, the boundary condition at x = 0 is not correctly set.

Your boundary condition is

D*dc/dx = h*(c-c_atm)

"pdepe" assumes a boundary condition of the form

pl + f*ql = 0

Rewriting your boundary condition in the required form gives

-h*(c-c_atm) + D*dc/dx = 0

, thus

pl = -h*(c-c_atm), ql = 1

Note the minus sign for pl.

You write that emission rates in the literature are related to area, not to volume. What surface area is meant here ? Can you give a link to a publication you are referring to ?

If you use google AI with the question "computing area related Root Respiration Emission Rate to volume related", you'll get a suggestion for conversion by a factor of 2/r where a root with cylindrical shape and radius r is assumed. Thus the S from the literature in mu-mol/m^2/s should be multiplied by epsilon * 2/r * 1e-6 to give unit mol/m^3/s where epsilon = volume root / (volume you use in the computation of c (usually volume of air)).

[Alister]

So in my model thus far it works less as a surface area but simply an exponential decay factor that scales the emissions by depth, according a study that defined the decay rate (Beta) for different types of plants. It's not a realistic representation as it doesn't even model a root system but the idea is to have some kind of source point represented in my model for the emissions. I will be expanding on this as I develop the model further so it more accurately estimates a root distribution.

I can certainly give an example publication Maize and Praire root CO2 Emissions.

Hope this helps.

[Torsten]

You wrote that emission rates in the literature are usually given in μmol m−2 s−1, not in mol/(m^3*s) as it is necessary in your code. Thus you have to convert a surface source (μmol m−2 s−1) to a volume source (mol/(m^3*s)). This has nothing to do with the exponential decay factor that scales the emissions by depth.