What about the sign of U1?

If I understand coreectly, at the interface you have

D1*( C1(i) - C1(i-1) ) - D2*( C2(i+1) - C2(i) ) = 0

where you have to set C1(i)=C2(i), isn't it?

Two things you may try to include you code to kill the interface instability:

Make the process of matching C1(i) and C2(i) iterative not master and slave, set a value for convergence. Consider Dirichlet ghost cell for the diffusion domain but its value come from AD domain, put both values in the iterative method so you get nice and smooth answer. Mesh refinement does not solve your problem here

Filippo Maria Denaro , thanks for you comments. U1 is advection velocity for the 1st layer. And yes, the second type BC at the interface with C1=C2.

the sign of U1 is negative?

Filippo Maria Denaro Yes, as written in the Fig. 1, it's more clear than in the text

Sorry but I haven’t found any infos about the value of U1

Kaveh Zamani thank you for your comments. I will try your suggestions. In fact, I used Crank-Nicolson scheme for this probleme. After discretization, I can get A*C(k+1) = C(k), k is time step, and A is the tridiagonal coefficient matrix. C1(i) and C2(i) in the vector C is the same thing C(i) and the BC in the interface is rewritten in A at the ith line. I can get the solution with the Thomas algorithm.

Filippo Maria Denaro Sorry, I misunderstood you. The values of the parameters are U1 = 2.e-4 cm/s; D1 = 2.e-4 cm2/s; D2 = 0.328 cm2/s; U3 = 0.024 cm/s. You can see D2>>D1, so the the gradient of C near the interface is very large. Thank you

Ok, if U1 is positive that makes sense of the Dirichlet BC.s at the lower boundary. Are you using a first order upwind for the convection?

Filippo Maria Denaro I used the space centre scheme for the convection, is it not good for this case?

Using central FD scheme (I suppose second order accurate) for both convection and diffusion terms you can have problems in the lack of monotonicity of the solution (that is you get numerical oscillations). To recover this problem you need to ensure that the Peclet cell number does not exceed the value of 2, this can be obtained by proper grid refinement.

Differently, a first order upwind ensures monotonicity.

Filippo Maria Denaro Actually, I havn't the problem of numerical oscillations, the Peclet number is 0.48 (=U1*dt/dz). My problem is the concentration (C) values of numerical solution (FD) is always smaller than analytical ones near the interface. I guess this is due to the great gradient of C near the interface. I have tried the upwind scheme, however, the error increases in contrast with the center scheme. I will try the higher accurate scheme near the interface with i-2, i-1, i, i+1 and i+2. Thanks a lot.

you are wrong, U1*dt/dz is not the Peclet number but the Courant number.

The cell Peclet number is U1*dz/D

Filippo Maria Denaro Sorry, I made a mistake. For U1*dz/D, in this case U1=D=2.e-4, so Peclet here is equal to dz, it's very small about 0.02. For the interface, we can write as -D1*C(i-1) + (D1+D2)*C(i) -D2*C(i+1) = 0, which is at the ith line of the matrix for crank-nicolson scheme or implicit scheme. I have tried increasing a little, e.g. 1/100 D1 for the first term and it becomes -1.01*D1*C(i-1) + (D1+D2)*C(i) -D2*C(i+1) = 0, the results become much better for all the column. 1/100 is not exact, but this shows maybe the problem is just the scheme for the interface.

At the interface are distinguishing C1 from C2?

No, actually, C1 and C2 are just for the mathematical model, for the numerical model, there is no C1 and C2, just C. After discretization, we can get the

A*C(k+1) = B*C(k), k is the time step. For the interface, there is only C(i). I have tried 5 points scheme instead of 3 points for the interface, the results improved. I will try 5 points scheme for all points and see the results.