NDSolve fails while solving PDE with Dancwerts BC

Question

The following code runs with BC1 (Dirichlet Type). However, it fails for Dancwerts type BC1.

Remove["Global`*"]
\[Rho]b = 651.52; \[Epsilon] = 0.59; Dia = 0.025; u = (.5*10^-6/60)/(Pi/4*Dia^2); L = 0.15; 
kc = 0.3217*10^-4; DL = 4.0259*10^-6; KL = .078; qm = 9*10^-3; C0 = 100*10^-6/10^-3;

Eq1 = D[Conc[t, z], t] + u* D[Conc[t, z], z] - DL* D[Conc[t, z], z, z] + (\[Rho]b/\[Epsilon])D[q[t, z], t] == 0;
Eq2 = D[q[t, z], t] == kc*(Conc[t, z] - q[t, z]/(KL (qm - q[t, z])));
IC1 = Conc[0, z] == 0;
IC2 = q[0, z] == 0;
(*BC1=(Conc[t,0]-DL/u *D[Conc[t,z],z])/.z\[Rule]0.00*) (*DANCKWERTS BC1*)
BC1 = Conc[t, 0];
BC2 = D[Conc[t, z], z] /. z -> L;

Soln = NDSolve[{Eq1, Eq2, IC1, IC2, BC1 == C0, BC2 == 0}, {Conc, q}, {t, 0, 10000}, {z, 0, L},Method -> {"Shooting", "StartingInitialConditions" -> {BC1 == C0}}, Method -> {StiffnessSwitching, Method -> {ExplicitRungeKutta, Automatic}}] 

Plot3D[Evaluate[Conc[t, z] /. Soln], {t, 0, 10000}, {z, 0.00, L}, PlotRange -> All] 
Plot[Evaluate[Conc[t, L] /. Soln], {t, 0, 10000}, PlotRange -> {0, C0}] 
Table[Flatten[{t, Evaluate[Conc[t, L] /. Soln]}], {t, 0, 10000, 500}] // TableForm

Answer 1

This problem can be solved with using Mathematica FEM. Let put u = -(.5*10^-6/60)/(Pi/4*Dia^2); and define Eq1 in accordance with Mathematica FEM as follows

Eq1 = D[Conc[t, z], t] - u*D[Conc[t, z], z] - 
   DL*D[Conc[t, z], z, z] + (\[Rho]b/\[Epsilon]) D[q[t, z], t];

Then Dancwerts type BC1 is

BC1 = (-u Conc[t, z] - DL*D[Conc[t, z], z]) /. z -> 0.00 

BC1 in FEM notation equals to BC1 = NeumannValue[-u (C0 - Conc[t, z]), z == 0]. From the other side BC2 = D[Conc[t, z], z] /. z -> L equals to BC2 = NeumannValue[0, z == L]. Therefore, code in FEM notation can be written as follows

Clear["Global`*"]
\[Rho]b = 651.52; \[Epsilon] = 0.59; Dia = 0.025; u = -(.5*10^-6/
      60)/(Pi/4*Dia^2); L = 0.15;
kc = 0.3217*10^-4; DL = 4.0259*10^-6; KL = .078; qm = 9*10^-3; C0 = 
 100*10^-6/10^-3;

Eq1 = D[Conc[t, z], t] - u*D[Conc[t, z], z] - 
   DL*D[Conc[t, z], z, z] + (\[Rho]b/\[Epsilon]) D[q[t, z], t];
Eq2 = D[q[t, z], t] == kc*(Conc[t, z] - q[t, z]/(KL (qm - q[t, z])));
IC1 = Conc[0, z] == 0;
IC2 = q[0, z] == 0;
(*BC1=Conc[t,0];*)
BC1 = NeumannValue[-u (C0 - Conc[t, z]), 
  z == 0];(*DANCKWERTS \
BC1=(Conc[t,0]-DL/u*D[Conc[t,z],z])/.z\[Rule]0.00*)

Soln = NDSolve[{Eq1 == BC1, Eq2, IC1, IC2}, {Conc, q}, {t, 0, 
   10000}, {z, 0, L}];

Visualization

{Plot3D[Evaluate[Conc[t, z] /. Soln], {t, 0, 10000}, {z, 0.00, L}, 
  PlotRange -> All, AxesLabel -> Automatic, PlotTheme -> "Marketing", 
  ColorFunction -> "Rainbow", MeshStyle -> White],
 Plot[Evaluate[{Conc[t, L], Conc[t, 0]} /. Soln], {t, 0, 10000}, 
  PlotRange -> {0, C0}, AxesLabel -> Automatic], 
 Plot[Evaluate[Table[Conc[t, z] /. Soln, {t, 1000, 10000, 1000}]], {z,
    0, L}, PlotRange -> {0, C0}, AxesLabel -> Automatic],
 Table[Flatten[{t, Evaluate[Conc[t, 0] /. Soln]}], {t, 0, 10000, 
    500}] // TableForm}

Note, that this solution computed on the mesh with 20 elements only,

Conc["ElementMesh"] /. Soln

Out[]= {NDSolve`FEM`ElementMesh[{{0., 
    0.15}}, {NDSolve`FEM`LineElement["<" 20 ">"]}]}  

Boundary conditions at z=0and z=L

{Show[Plot[(-u Conc[t, 0] - DL Derivative[0, 1][Conc][t, 0]) /. 
    Soln, {t, 0, 10000}, PlotRange -> All], 
  ListPlot[Table[{t, -u C0}, {t, 1000, 10000, 1000}], 
   PlotStyle -> Red]], 
 Plot[Derivative[0, 1][Conc][t, L] /. Soln, {t, 0, 10000}, 
  PlotRange -> All]}

We can also compute solution with DirichletCondition[] as follows

BC1 = Conc[t, 0];
BC2 = D[Conc[t, z], z] /. z -> L;

Sol0 = NDSolve[{Eq1 == 0, Eq2, IC1, IC2, 
    DirichletCondition[Conc[t, z] == C0, z == 0]}, {Conc, q}, {t, 0, 
    10000}, {z, 0, L}];


{Plot3D[Evaluate[Conc[t, z] /. Sol0], {t, 0, 10000}, {z, 0.00, L}, 
  PlotRange -> All],
 Plot[Evaluate[{Conc[t, L], Conc[t, 0]} /. Sol0], {t, 0, 10000}, 
  PlotRange -> {0, C0}, AxesLabel -> Automatic], 
 Plot[Derivative[0, 1][Conc][t, L] /. Sol0, {t, 0, 10000}],
 Table[Flatten[{t, Evaluate[Conc[t, L] /. Sol0]}], {t, 0, 10000, 
    500}] // TableForm}