Mass Flux Output Issue for a 1D Diffusion and Linear Sorption Problem

[Girish Kumar]

Hello Glenn,

Hope all is well! I'm running the attached simple 1D diffusion problem with linear sorption isotherm in PFLOTRAN. The model is a 1D column of clay that is 8 cm in length, porosity of 0.5, rock density (soil particle density) of 2650 kg/m3. The top and bottom boundary are dirichlet boundary conditions with the bottom at a known fixed concentration of 9.13 mM and the top is at a zero concentration for the species in question. I assign a Kd value of 500 in the default units (kgw/m3_bulk) for the species to the clay. This corresponds to a retardation factor of 2. 

When I run the simulation for a certain time, print the instantaneous flux at the top boundary and compare it with an analytical solution for this problem with the same parameters, they do not match. However, when I divide the instantaneous flux produced by PFLOTRAN by the retardation factor, it matches with the analytical solution. So, I wonder if there's any error in the flux output or if I'm missing something? Since the retardation affects the concentration profile and that concentration affects the gradient at the top and thereby the flux at that boundary, I'd expect the effect of retardation to also translate to the flux. But that doesn't seem to be happening.

The analytical solution is equation 4.24 from the Crank book as shown below and it is the solution for cumulative flux (mol/m^2). I take its time derivative which gives me instantaneous mass flux in mol/m^2/s:

To incorporate retardation into the analytical solution I just divide the diffusion coefficient by the retardation factor. This assumption is valid for linear equilibrium sorption (Kd), with a retardation factor associated with only the accumulation term in the governing equation. Note that I multiply the instantaneous flux by porosity to make this solution applicable to porous media.

I do not expect you to create a script for the analytical solution, but if you can check that there's no error in the mass flux output at the boundaries that'd be helpful. 

I have attached the input file and the database file (hanford.dat) for your reference. 

Thanks,

Girish

--

E-mail: girish...@gmail.com

[Hammond, Glenn E]

Girish,

Could you try using a retardation factor other than 2 and let me know how your results compare to the analytical solution?

Our boundary fluxes are calculated as first order, meaning the gradient is evaluated over half the cell width. In contrast, internal fluxes use a full cell width—half on each side—assuming uniform resolution. This half-distance approach may be causing the boundary flux to double.

To further investigate, you can add 2ND_ORDER_BOUNDARY_CONDITION in the GRID block. 

https://documentation.pflotran.org/user_guide/cards/subsurface/grid_card.html?highlight=2nd_order_boundary_condition

This adjustment sets the boundary connection spacing to the full cell width by essentially introducing a fictitious node half a cell width outside the domain. I am wondering what you observe with this setting.

Regards,

Glenn

From: pflotra...@googlegroups.com <pflotra...@googlegroups.com> on behalf of Girish Kumar <girish...@gmail.com>
Date: Thursday, April 23, 2026 at 5:49 PM
To: pflotra...@googlegroups.com <pflotra...@googlegroups.com>
Subject: [pflotran-users: 8775] Mass Flux Output Issue for a 1D Diffusion and Linear Sorption Problem

Check twice before you click! This email originated from outside PNNL.

--
You received this message because you are subscribed to the Google Groups "pflotran-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pflotran-user...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/pflotran-users/CADZRAywqn%3D_7UQbPEtDkzJ%3DFULZX5td9G5O26-H0ouN_grNgVg%40mail.gmail.com.

[Girish Kumar]

Hi Glenn,

I see the nuance. Since the analytical solution directly provides the equation without the need for "dx" in dC/dx, I wrote a finite difference code for this problem (I verified that the FDM solution matches the analytical solution) with a discretization such that dx at the boundary is same as the dx in my PFLOTRAN file. The issue still persists for any other retardation factor. 

I even tried the 2nd order boundary condition and the issue remains. Please advise.

Below is the FDM code in Matlab (I don't know how I can attach files here while replying).

%%%%%%MATLAB Code for FDM Solution%%%%%%%%%%%%%%%

clc; clear;

%% PARAMETERS

L = 0.08; % domain length [m]

Nx = 80; % number of cells

dx = L / Nx; % grid spacing = 1 mm per cell

x = linspace(0, L, Nx+1)'; % 81 nodes

poro = 0.5; % porosity

R = 2.0; % Retardation factor

tort = 1.0; % tortuosity (right now not considering it)

D = tort*1.0e-9/R; % diffusion coefficient [m2/s]

C_left = 1200.0; % ppm = mg/L (corresponds to 9.13 mM for TCE)

C_right = 0; % ppm = mg/L

A_cm2 = 9.62; % cross-sectional area [cm2]

A = A_cm2 * 1e-4; % convert to m2

t_final_days = 24; % simulation time

t_final = t_final_days * 24*3600; % seconds

%% STABILITY CONDITION

dt_max = dx^2 / (2*D);

dt = 0.5 * dt_max; % a safe dt

Nt = ceil(t_final / dt);

fprintf("dx = %g m\n", dx);

fprintf("Stable dt <= %.3e s, using dt = %.3e s\n", dt_max, dt);

fprintf("Total time steps = %d\n", Nt);

%% INITIAL CONDITION

C = zeros(Nx+1,1); % ppm = mg/L

%% STORAGE FOR FLUX

flux_mg_per_day = zeros(Nt,1);

time_days = zeros(Nt,1);

C_at_node80 = zeros(Nt,1); % storage array for node 80 concentration

%% TIME STEPPING

for n = 1:Nt

C_old = C;

% Boundary conditions

C_old(1) = C_left;

C_old(end) = C_right;

% Update interior nodes

for i = 2:Nx

C(i) = C_old(i) + D*dt/dx^2 * (C_old(i+1) - 2*C_old(i) + C_old(i-1));

end

% Reapply BCs

C(1) = C_left;

C(end) = C_right;

%% ---- INSTANTANEOUS MASS FLUX AT x = L ----

% Concentration gradient (mg/L per m)

dCdx = (C(end) - C(end-1)) / dx; % backward difference

% Fick's law flux (mg/L * m/s)

J = -D * poro * dCdx*1000; % mg/m^3 * m/s = mg/m^2/s

% Multiply by area (m2) → mg/s

mass_flux_mg_s = J * A;

% Convert to mg/day

mass_flux_mg_day = mass_flux_mg_s * 86400;

% Store

flux_mg_per_day(n) = mass_flux_mg_day;

time_days(n) = n * dt / 86400;

%% ---- STORE CONCENTRATION AT NODE 80 ----

C_at_node80(n) = C(Nx); % node just left of boundary

end

%% FINAL CONCENTRATION PLOT

figure;

plot(x*100, C, 'b-', 'LineWidth', 2);

xlabel('Distance (cm)');

ylabel('Br^- concentration (ppm)');

title('1D Diffusion After 23 Days (Numerical FD, 80 cells)');

grid on;

%% FLUX PLOT

figure;

set(gcf,'WindowState','maximized')

clf

plot(time_days, flux_mg_per_day./0.017, 'r-', 'LineWidth', 2);

xlabel('Time (days)');

ylabel('Instantaneous flux (mg/day)');

title('Mass Flux at x = L (outlet face)');

grid on;

 

Thanks,

Girish

[Hammond, Glenn E]

Girish,

Thank you for conducting additional simulations to rule out boundary discretization as the source of the issue.

At this stage, I suspect the problem is more conceptual than mathematical, but I want to confirm this for myself.

My current understanding is that, according to the equation (above), dividing by the retardation factor should reduce both flow velocity and diffusion by half. However, for a given aqueous concentration, both internal and boundary fluxes remain unaffected by sorption, since advective and dispersive fluxes depend solely on the aqueous phase concentration and concentration gradient, respectively. Meanwhile, half the mass is retained in the sorbed phase. This leads me to wonder if it is actually the concentration used in the flux calculation in the equation above that should  be halved.

Thought process: From a mass perspective, for R=2 each grid cell divides its stored mass evenly between the aqueous and sorbed phases. Only the mass in the aqueous phase is capable of undergoing advection and dispersion. If you were to add an INTEGRAL_FLUX at the last internal connection (other side of the boundary cell), there's a strong possibility that the resulting flux would also appear doubled….

Would you agree that dividing by the retardation factor, R, causes the concentration C in the accumulation term (the time derivative) to reflect the total mass present within the cell? Consequently, should the aqueous concentration used in the flux calculation be half of the concentration accounted for in the accumulation term?

I will also take some time to consider this further, and welcome input from anyone else who would like to share their perspective.

Glenn

To view this discussion visit https://groups.google.com/d/msgid/pflotran-users/f40bc696-ab6e-4a8a-96bc-b04c102274c7n%40googlegroups.com.

[Girish Kumar]

Hello Glenn,

Thanks for taking the time to mull on this. Based on my understanding, retardation does not alter how flux is computed, i.e., flux would still be J = D*dC/dx, and not D/R*dC/dx. Instead, retardation appears in the accumulation term, slowing the evolution of Caq, i.e, dC/dt is effectively reduced by a factor of R. Because of this, I think the Caq used in the flux calculation should not be halved (for R=2 or otherwise). Caq evolves according to the retarded transport equation rather than being rescaled by R. Any differences in flux between retarded and non-retarded cases arise indirectly through differences in the concentration field.

So, I think I was making a mistake (and in fact double-counting retardation) by introducing 1/R for the flux term in my analytical and FDM solution, and thereby underestimating the aqueous concentrations. I think the flux printed by PFLOTRAN are actually correct. After digging through the source code, I see that PFLOTRAN computes boundary flux using the concentration values resulting from retardation near the boundary and multiplies with effective diffusivity, which I think is correct. 

Let me know if this solidifies your conceptual understanding and that we're on the same page.

Regards,

Girish