Update: Mass Transport Tutorial

In my answer, I adapted the Heat Transfer Tutorial to handle Mass Transport. There is a new and fairly extensive Mass Transport Tutorial that you can use directly. In my answer 227821, I extended the Mass Transport functions in the tutorial to include axisymmetric and optional porosity.

Update: Extended Answer

To use the FEM method it is good to cast your equations into coefficient form as shown FEM Tutorial.

$$frac{{{partial ^2}}}{{partial {t^2}}}u + dfrac{partial }{{partial t}}u + nabla cdotleft( { - cnabla u - alpha u + gamma } right) + beta cdotnabla u + au - f = 0$$

The benefits of coefficient form include:

• Standardization so others can quickly interpret your equations
• Avoidance of sign errors
• Nice one-to-one mapping with other FEM codes such as COMSOL

For advective-diffusive problems with and without source terms, the Heat transfer tutorial is a good place to start to learn how to set up these types of problems. It shows how to properly setup the PDE operators in coefficient form. The TimeHeatTransferModel function in the Heat Transfer Tutorial can help guide you on how to setup the problem.

Create Fluid and Solid PDE Operators Using Heat Transfer Tutorial

The following shows how one can setup your coupled system of interphase transport of the fluid and solid using the TimeHeatTransferModel :

(* From Vitaliy Kaurov for nice display of operators *)
pdConv[f_] := 
 TraditionalForm[
  f /. Derivative[inds__][g_][vars__] :> 
    Apply[Defer[D[g[vars], ##]] &, 
     Transpose[{{vars}, {inds}}] /. {{var_, 0} :> 
        Sequence[], {var_, 1} :> {var}}]]
(*PDEModels/tutorial/HeatTransfer/HeatTransferVerificationTests#
463435833*)
ClearAll[HeatTransferModel]
HeatTransferModel[T_, X_List, k_, ρ_, Cp_, Velocity_, Source_] :=
  Module[{V, Q, a = k},
  V = If[Velocity === "NoFlow", 
    0, ρ*Cp*Velocity.Inactive[Grad][T, X]];
  Q = If[Source === "NoSource", 0, Source];
  If[ FreeQ[a, _?VectorQ], a = a*IdentityMatrix[Length[X]]];
  If[ VectorQ[a], a = DiagonalMatrix[a]];
  (* Note the - sign in the operator *)
  a = PiecewiseExpand[Piecewise[{{-a, True}}]];
  Inactive[Div][a.Inactive[Grad][T, X], X] + V - Q]
TimeHeatTransferModel[T_, TimeVar_, X_List, k_, ρ_, Cp_, 
  Velocity_, Source_] := ρ*Cp*D[T, {TimeVar, 1}] + 
  HeatTransferModel[T, X, k, ρ, Cp, Velocity, Source]
(* Create Parametric PDE Operators for Fluid and Solid *)
(* Include Small Diffusive Term to Fluid Otherwise MMA Might Complain 
about Being Convection Dominated *)
parmfop = 
  TimeHeatTransferModel[c[t, x], t, {x}, {d}, 1, 
   1, {u}, -(ρ/ϵ) Q];
parmsop = 
  TimeHeatTransferModel[q[t, x], t, {x}, {0}, 1, 1, "NoFlow", Q];
parmfop // pdConv
parmsop // pdConv

Slight rearrangement shows the operators follow coefficient form:

$$begin{matrix} mfrac{{{partial ^2}}}{{partial {t^2}}}u & + dfrac{partial }{{partial t}}u & + nabla cdotleft( { - cnabla u - alpha u + gamma } right) u & + beta cdotnabla u & + au & - f & =0 & (1) & & & & & & & + frac{{partial c(t,x)}}{{partial t}} & + {nabla _{left{ x right}}}cdotleft( {left( {begin{array}{*{20}{c}} { - d} end{array}} right).{nabla _{left{ x right}}}cleft( {t,x} right)} right) & + left{ u right}.{nabla _{left{ x right}}}cleft( {t,x} right) & & + frac{{rho Q}}{varepsilon } & =0 & (2) & & & & & & & + frac{{partial q(t,x)}}{{partial t}} & & & & - Q & =0 & (3) end{matrix}$$

New Workflow

The following workflow uses the PDE operators defined above and produces an annotated plots that can be scrolled frame-by-frame with manipulate.

(* Define Parameters *)
eps = 0.373;
ueps = 6.66;
c0 = 892;
Ks = 0.0872;
rho = 390.04;
qm = 7.0136;
b = 0.05313;
tend = 60;
xend = 22;
Qsource = -Ks*(q[t, x] - (qm*b*c[t, x])/(1 + b*c[t, x]));
timeinc = 0.25; (* Output times *)
(* Critical Breakthrough Fraction *)
btfrac = 0.05; 
(* Use CombinePlots For Secondary Axis *)
cp = ResourceFunction["CombinePlots", "Function"];
(* Setup PDE System *)
pdef = (parmfop == 0) /. {ϵ -> eps, u -> ueps, ρ -> rho,
     d -> 1, Q -> Qsource};
pdes = (parmsop == 0) /. {Q -> Qsource};
dc = DirichletCondition[c[t, x] == c0 (1 - Exp[-40 t]), x == 0];
icf = c[0, x] == 0;
ics = q[0, x] == 0;
(* Solve PDE *)
{cifn, qifn} = 
  NDSolveValue[{pdef, pdes, dc, icf, ics}, {c, q}, {t, 0, tend}, {x, 
    0, xend}, 
   Method -> {"MethodOfLines", "TemporalVariable" -> t, 
     "SpatialDiscretization" -> {"FiniteElement", 
       "MeshOptions" -> {"MaxCellMeasure" -> 0.5}}}];
(* Plot Solutions *)
(* Define Plot Function Using Combine Plots *)
bfn = Function[{t}, 
   If[cifn[t, xend] >= btfrac c0, xend, 
    x /. FindRoot[cifn[t + 0.01, x] - btfrac c0, {x, 0, xend}, 
      Method -> "Brent"]]];
bttemp = With[{xbt = #}, 
    If[xbt == xend, "Sat", 
     StringTemplate["x=``"][ToString@PaddedForm[xbt, {4, 2}]]]] &;
pltfn = Function[{t}, cp[
    With[{xbt = bfn[t]}, 
     Plot[Callout[cifn[t, x], bttemp[xbt], xbt, RoundingRadius -> 5, 
       Appearance -> "SlantedLabel"], {x, 0, xend}, 
      PlotRange -> {{-2, xend + 1}, {-10, 1.1 c0}}, Frame -> True, 
      FrameLabel -> {"x", "c"}, 
      PlotLabel -> 
       Style[StringTemplate[
          "Fluid and Solid Concentrations @ time=``"][
         ToString@PaddedForm[t, {4, 2}]], 12], 
      Epilog -> {Green, Thick, Dashed, 
        Line[{{xbt, -10}, {xbt, 1.1 c0}}]}]],
    Plot[qifn[t, x], {x, 0, xend}, PlotPoints -> 100, 
     PlotRange -> {0, 7},
     Frame -> True, FrameStyle -> Red, PlotStyle -> Red, 
     FrameLabel -> {"x", "q"}],
    "AxesSides" -> "TwoY"
    ], Listable];
(* Create Image List and Manipulate *)
imgs = Rasterize[#, ImageSize -> 340] & /@ 
   pltfn[Range[0, tend, timeinc]];
Manipulate[imgs[[i]], {i, 1, tend/timeinc, 1}, 
 ControlPlacement -> Top]

Original Answer

I would cast the problem using the following workflow:

(* Define Parameters *)
eps = 0.373;
ueps = 6.66;
c0 = 892;
Ks = 0.0872;
rho = 390.04;
qm = 7.0136;
b = 0.05313;
tend = 60;
xend = 22;
Qsource = -Ks*(q[t, x] - (qm*b*c[t, x])/(1 + b*c[t, x]));
(* Use CombinePlots For Secondary Axis *)
cp = ResourceFunction["CombinePlots"];
(* Define Parametric PDE operators *)
(* Include Small Diffusive Term to Fluid Otherwise MMA Might Complain 
about Being Convection Dominated *)
parmfop = 
  D[c[t, x], t] +  D[-d D[c[t, x], x], x] + u  D[c[t, x], x] + (
   Q ρ)/ϵ;
parmsop = D[q[t, x], t] - Q;
(* Setup PDE System *)
pdef = (parmfop == 0) /. {ϵ -> eps, u -> ueps, ρ -> rho,
     d -> 1, Q -> Qsource};
pdes = (parmsop == 0) /. {Q -> Qsource};
dc = DirichletCondition[c[t, x] == c0 (1 - Exp[-40 t]), x == 0];
icf = c[0, x] == 0;
ics = q[0, x] == 0;
(* Solve PDE *)
{cifn, qifn} = 
  NDSolveValue[{pdef, pdes, dc, icf, ics}, {c, q}, {t, 0, tend}, {x, 
    0, xend}, 
   Method -> {"MethodOfLines", "TemporalVariable" -> t, 
     "SpatialDiscretization" -> {"FiniteElement", 
       "MeshOptions" -> {"MaxCellMeasure" -> 0.5}}}];
(* Plot Solutions *)
(* Define Plot Function Using Combine Plots *)
pltfn = Function[{t}, cp[
    Plot[cifn[t, x], {x, 0, xend}, PlotRange -> {0, 1.1 c0}, 
     Frame -> True, FrameLabel -> {"x", "c"}],
    Plot[qifn[t, x], {x, 0, xend}, PlotPoints -> 100, 
     PlotRange -> {0, 7},
     Frame -> True, FrameStyle -> Red, PlotStyle -> Red, 
     FrameLabel -> {"x", "q"}],
    "AxesSides" -> "TwoY"
    ], Listable];
(* Create Image List and Animate *)
imgs = Rasterize /@ pltfn[Range[0, tend, 0.25]];
ListAnimate[imgs, ControlPlacement -> Top]