*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
```