# Mineral dissolution and solute transport around a solid

I am trying to simulate solute transport of acid (HCl) and consequent mineral dissolution around a grain (calcite).

The governing equation for transport is the advection-diffusion equation, given as:

where C is the concentration, v is the flow velocity, D is the diffusion coefficient. At a constant flow velocity, the above is discretized using an FTCS formulation as:

The dissolution occur at the fluid/solid interface – and can be modelled as a Robin type Bc, which is described as:

where r is the mineral rate of dissolution, ε is the stoichiometric coefficient of species in the dissolution process (this is equal to 1 in the reaction considered), k_(H+) is the rate constant, γ_(H+) is the activity coefficient, and c_(H+) is the molar concentration of reacting specie (in this case Hydrogen or chloride ion).

I have discretized Eq. (4) as:

This is my code so far:

clc; clear; close all;

%%%%%%%%%% Specify inputs
CircleDiam = 40;        %in pixels
Co    = 0;          %initial conc. in the domain [mol/L]
Cin   = 0.01;       %conc. of injected fluid [mol/L]

Lx    = 0.1/100;        %Length of domain [m]
nx    = 200;        %spatial gridpoints in x
dx    = Lx/(nx-1);      %Length step size [m]
Ly    = 0.05/100;       %Width of domain [m]
ny    = 100;        %spatial gridpoints in y
dy    = Ly/(ny-1);      %Length step size [m]

T     = 5/(3600*24);    %Simulation time [days]
nt    = 8000;           %shifts
dt    = T/nt;       %Time step [days]

% Flow
u = 103.68;         %Velocity in x direction [m/day]
v = 0;              %Velocity in y direction [m/day]
De = 8.64e-04;      %Dispersion coeff. [m2/day]
betaX  = u*dt/dx;
betaY  = v*dt/dy;
gammaX = De*dt/(dx^2);
gammaY = De*dt/(dy^2);

%%%%%%%%%% Create image with solid object

% obtain full output grids from grid vectors
[Colgrids, Rowgrids] = meshgrid(1:nx, 1:ny);

% create a logical mask for the circle by specifying the center and diameter of the circle.
centerX = 0.5 + (nx/2);
centerY = 0.5 + (ny/2);

% obtain image from: ( (y-y0)^2 + (x-x0)^2 )  <= r^2, where (y0,x0) is the centre point of circle
SolidImg = (Rowgrids - centerY).^2 + (Colgrids - centerX).^2 <= radius.^2;

% change from logical to numeric labels. Also, transpose matrix to conform with the conc. matrix
P   = double(SolidImg');
%figure, imshow(~SolidImg, [], 'InitialMagnification','fit'); box on;

%%%%%%%%%% Find the solid boundaries
BW = im2bw(P);
BW_filled  = imfill(BW,'holes');
boundaries = bwboundaries(BW_filled);

for k = 1:size(boundaries, 1)
bdr_matrix = boundaries{k};
end

row = bdr_matrix(:,1);
col = bdr_matrix(:,2);

% Gridblocks
x = 0:dx:Lx;
y = 0:dy:Ly;
t = 0:dt:T;
[X,Y] = meshgrid(x, y);

% specify initial conditions
C        = zeros(nx, ny, nt+1);
C(:,:,1) =  Co;              %Initial condition

% specify Reaction parameters
kh     = 77004.081;     % [mol/m2/day]
gammaH = 1;             % [L/mol]

%iterate finite difference equations
for k = 1:nt
for j = 2:ny-1
for i = 2:nx-1
if P(j,i)==1    %Solid pixels
C(j,i,k+1) = C(j,i,k);

% Apply the robin bc at the fluid/solid interface
if ismember(P(j,i), P(row, col))
C(j,i,k+1) = C(j,i,k)...
- 1000*0.5*De*(C(j,i+1,k)-C(j,i-1,k)/dx)/(kh*gammaH)...
- 1000*0.5*De*(C(j+1,i,k)-C(j-1,i,k)/dy)/(kh*gammaH); % multiply by 1000L to make dimensionally consistent with mol/L
end
else
C(j,i,k+1) = C(j,i-1,k)*(betaX/2+gammaX) + C(j,i+1,k)*(gammaX-betaX/2)...
+ C(j-1,i,k)*(betaY/2+gammaY) + C(j+1,i,k)*(gammaY-betaY/2)...
+ C(j,i,k)*(1-2*gammaX-2*gammaY);
end
end
end

% Insert boundary conditions
C(1,:,k+1)   = C(2,:,k+1);      % south
C(end,:,k+1) = C(end-1,:,k+1);  % north
C(:,1,k+1)   = Cin;             % west
C(:,end,k+1) = C(:,end-1,k+1);  % east
end


My result (see image) does not look like what I would expect (see second image). Ideally, I expect to see something like the above image, and so I am wondering where I might be getting it wrong. I have assumed that both velocity and diffusion coefficient are zero at the solid cells, and I have also prescribed a robin type bc at the fluid/solid interface, however, I get the same outcome when I run the code without eq. (5). I want to believe that I am not doing this the right way. I am not completely sure if I have discretized Eq. (5) the right way, and I am wary of how I have implemented it (i.e., eq. 5) in the code. I thus would very much appreciate if someone can help me out here or point me in the right direction.

PS: The axes are the same in both image. The difference is that mine is in pixels while the other is in cm.

Thank you in anticipation.

EDIT

Thank you so much for your help @Sthavishtha Bhopalam . I am indeed very grateful.

1. I obtained the velocity field from LBM. It’s a pity I didn’t mention this initially. I have now included the appropriate lines of code (without the LBM though) as shown below - and I'd appreciate if you can confirm the correctness.

The following comes just before the loop:

 % calculate flow velocity field
[u,v]  = FlowVelocitySolver(nx,ny,1000);  % iter = 1000;


And then the values of u and v at every point in the grid is utilized in the advection term.

1. I have run the code with your amendments – and I got the expected shape as in Molins et al. (see the new figure), however:

• I see you changed the inlet velocity from 103.68 m/day (which is 0.12 cm/s) to 10.368 m/day (which is 0.012 cm/s), the former being the same with Molins et al.

• I noticed that numerical instability occurs at the fluid-solid interface when I use 103.68 m/day. This happens even after I adjusted the time step value to match the stability criteria. Is this the reason you changed to 10.368 m/day?

1. I was looking to use other kinetic formulations. For example, see pages 68-69 of: https://biblio.ugent.be/publication/8071415/file/8071428.pdf

Would it be right to discretize Eqs. 4.10, given by:

$$-D\textbf n.\nabla c = \frac{V}{A_{s}}k_{het}(c – c_{eq})$$

as:

Many thanks once again.

• I tried running this in Octave, an open source Matlab clone, but nothing came out. And I don't want to convert it to Python by hand to poke at it. I'd like to help, but I don't want to buy Matlab. Your choice of programming language has made it difficult for me to help you :-( Dec 2, 2021 at 7:16
• @Richard, my bad I am still a python learner at the moment. Thank you for even trying to help. Cheers! Dec 2, 2021 at 11:05

Firstly, there are some mistakes in your discretized formulation and the approach you employ to solve this problem:

• The flow velocity in the advection diffusion Eq. (1) is to be computed from the fluid continuity (COM) and momentum equations (COLM). In other words, you cannot assume the flow velocity to be constant in the entire domain, i.e., you employ constant variables u and v in the code. The correct result (which you extract from ref. [1]) is intuitively identical to the flow over a two-dimensional cylinder at a low Reynolds number. So, this necessitates the need to solve COM and COLM.
• To ensure numerical stability, the advective part of Eq. (1) has to be discretized by an upwind scheme. So, the discrete form would be (for the purpose of illustration here, I assume $$u > 0, v > 0$$)

$$\frac{c_{i,j}^{n+1} - c_{i,j}^{n}}{\Delta t} = D \left(\frac{c_{i+1,j}^n - 2c_{i,j}^n + c_{i-1,j}^n}{(\Delta x)^2} + \frac{c_{i,j+1}^n - 2c_{i,j}^n + c_{i,j-1}^n}{(\Delta y)^2}\right) - \frac{(uC)_{i,j}^n - (uC)_{i-1,j}^n}{\Delta x} - \frac{(vC)_{i,j}^n - (vC)_{i,j-1}^n}{\Delta y}$$

• After obtaining the discrete equation, you must appropriately choose the input parameters to ensure the satisfaction of the CFL criterion, i.e., positivity of the coefficient of $$c_{i,j}^n$$.

• The discrete form Eq. (5) of the Robin BC Eq. (4) needs to be changed. The LHS of the discrete form would be (I just show forward difference here, but you could also use central differencing)

$$-D\left(n_x \frac{C_{i,j}^n - C_{i-1,j}^n}{\Delta x} + n_y \frac{C_{i,j}^n - C_{i,j-1}^n}{\Delta y}\right) = RHS$$

where $$n_x, n_y$$ are the unit normal vectors to the fluid-solid interface. You may have to double check the direction of these vectors from appropriate references. Also, I think that the term $$c_{H+}$$ in the RHS is the concentration at the fluid-solid interface, which would not be equivalent to $$c_{i,j}^{n+1} - c_{i,j}^{n}$$.

• I have made some of the aforementioned changes to your script, except of solving for COM and COLM. As you may see, the solution still partially resembles your plots, primarily because of neglecting these equations.
1. Molins, Sergi; Soulaine, Cyprien; Prasianakis, Nikolaos I.; Abbasi, Aida; Poncet, Philippe; Ladd, Anthony J. C.; Starchenko, Vitalii; Roman, Sophie; Trebotich, David; Tchelepi, Hamdi A.; Steefel, Carl I., Simulation of mineral dissolution at the pore scale with evolving fluid-solid interfaces: review of approaches and benchmark problem set, ZBL07420582.
    clc; clear; close all;
% parameters slightly modified as per the paper
%%%%%%%%%% Specify inputs
CircleDiam = 40;    %in pixels
Co    = 0;          %initial conc. in the domain [mol/L]
Cin   = 0.01;       %conc. of injected fluid [mol/L]

Lx    = 0.1/100;    %Length of domain [m]
nx    = 200;        %spatial gridpoints in x
dx    = Lx/(nx-1);  %Length step size [m]
Ly    = 0.05/100;   %Width of domain [m]
ny    = 100;        %spatial gridpoints in y
dy    = Ly/(ny-1);  %Length step size [m]

T     = 500/(3600*24);%Simulation time [days]
nt    = 5E4;        %shifts
dt    = T/nt;       %Time step [days]

% Flow
% u = 103.68;       %Velocity in x direction [m/day]
u = 10.368;
v = 0;              %Velocity in y direction [m/day]
De = 8.64e-06;      %Dispersion coeff. [m2/day]
betaX  = u*dt/dx;
betaY  = v*dt/dy;
gammaX = De*dt/(dx^2);
gammaY = De*dt/(dy^2);

%%%%%%%%%% Create image with solid object

% obtain full output grids from grid vectors
[Colgrids, Rowgrids] = meshgrid(1:nx, 1:ny);

% create a logical mask for the circle by specifying the center and diameter of the circle.
centerX = 0.3 + (nx/2);
centerY = 0.3 + (ny/2);

% obtain image from: ( (y-y0)^2 + (x-x0)^2 )  <= r^2, where (y0,x0) is the centre point of circle
SolidImg = (Rowgrids - centerY).^2 + (Colgrids - centerX).^2 <= radius.^2;

% change from logical to numeric labels. Also, transpose matrix to conform with the conc. matrix
P   = double(SolidImg');
%figure, imshow(~SolidImg, [], 'InitialMagnification','fit'); box on;

%%%%%%%%%% Find the solid boundaries
BW = im2bw(P);
BW_filled  = imfill(BW,'holes');
boundaries = bwboundaries(BW_filled);

for k = 1:size(boundaries, 1)
bdr_matrix = boundaries{k};
end

% gxn and gyn have been accidentally swapped
gxn = gx;
gyn = gy;

for i = 1:nx
for j = 1:ny
gxn(i,j) = gx(i,j)/gmag(i,j);
gyn(i,j) = gy(i,j)/gmag(i,j);
end
end

row = bdr_matrix(:,1);
col = bdr_matrix(:,2);

% Gridblocks
x = 0:dx:Lx;
y = 0:dy:Ly;
t = 0:dt:T;
[X,Y] = meshgrid(x, y);

% specify initial conditions
C        = zeros(nx, ny);
C(:,:)   = Co;              %Initial condition
Cnew     = C;

% specify Reaction parameters
kh     = 77004.081;     % [mol/m2/day]
gammaH = 1;             % [L/mol]

%iterate finite difference equations
for k = 1:nt
for j = 2:ny-1
for i = 2:nx-1
if P(i,j)==1    %Solid pixels
Cnew(i,j) = 0.;

% Apply the robin bc at the fluid/solid interface
for l = 1:size(row,1)
if (i == row(l)) && (j == col(l))
temp = ...
- 1000*De*.5*gyn(i,j)*(C(i+1,j)-C(i-1,j)/dx)/(kh*gammaH)...
- 1000*De*.5*gxn(i,j)*(C(i,j+1)-C(i,j-1)/dy)/(kh*gammaH); % multiply by 1000L to make dimensionally consistent with mol/L
Cnew(i,j) = temp;
end
end
else
Cnew(i,j) = C(i-1,j)*(betaX/2+gammaX) + C(i+1,j)*(gammaX-betaX/2)...
+ C(i,j-1)*(betaY/2+gammaY) + C(i,j+1)*(gammaY-betaY/2)...
+ C(i,j)*(1-2*gammaX-2*gammaY);
end
end
end

% Insert boundary conditions
Cnew(:,1)   = Cnew(:,2);      % south
Cnew(:,end) = Cnew(:,end-1);  % north
Cnew(1,:)   = Cin;            % west
Cnew(end,:) = Cnew(end-1,:);  % east

k = k + 1;
C = Cnew;
end

• you should post the questions/comments related to my answer here, instead of editing in the original post. Reg. the flow velocity, have you tried changing the diffusion coefficient to $8.64e-5$, as used in the paper. you need to do the same for other parameters you employ in your lbm code. Dec 14, 2021 at 23:54
• also, I still spot some mistakes in the changes you made in the code. $(uC)_{i-1,j}$ is not equivalent to $u_{i,j}C_{i-1,j}$, as $u$ is varying (obtained from your flow solver). Dec 15, 2021 at 0:00
• reg. eqn. 4.10: it looks like $c - c_{eq}$ should actually be $c_{\Gamma} - c_{eq}$, where $c_{\Gamma}$ denotes the concentration at the fluid-solid interface. and $c_{\Gamma}$ should be somehow linked to the solid concentration; but i am not very sure. Dec 15, 2021 at 0:29
• Thanks again for your response. You're right, I shouldn't have added anything to your code. So, I have now removed my additions as appropriate. Furthermore, I didn't realize you changed De from Molins. In any case, I thought the original value should be 8.64e-4 ? Or I'm I missing something? Dec 15, 2021 at 13:01
• Thank you also for your comment regarding eqn. 4.10. As pertaining the expansion of $(uC)_{i-1,j}$, please what should be the correct thing? Dec 15, 2021 at 13:03