### 2.25.7.5. Crank-Nicholson Solution[¶](#crank-nicholson-solution "Permalink to this headline")

Equation [(2.25.1)](#equation-24-1) is solved using a Crank-Nicholson solution ([Press et al., 1992](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/References/CLM50_Tech_Note_References.html#pressetal1992)), which combines fully explicit and implicit representations of the mass balance. The fully explicit decomposition of equation [(2.25.1)](#equation-24-1) can be written as

(2.25.15)[¶](#equation-24-15 "Permalink to this equation")\\\[\\frac{R\_{j}^{n+1} C\_{j}^{n+1} -R\_{j}^{n} C\_{j}^{n} }{\\Delta t} =\\frac{1}{\\Delta x\_{j} } \\left\[\\frac{D\_{p1}^{n} }{\\Delta x\_{p1}^{} } \\left(C\_{j+1}^{n} -C\_{j}^{n} \\right)-\\frac{D\_{m1}^{n} }{\\Delta x\_{m1}^{} } \\left(C\_{j}^{n} -C\_{j-1}^{n} \\right)\\right\]+S\_{j}^{n} ,\\\]

where _j_ refers to the cell in the vertically discretized soil column (increasing downward), _n_ refers to the current time step, \\(\\Delta\\)_t_ is the time step (s), _p1_ is _j+½_, _m1_ is _j-½_, and \\(S\_{j}^{n}\\) is the net source at time step _n_ and position _j_, i.e., \\(S\_{j}^{n} =P\\left(j,n\\right)-E\\left(j,n\\right)-A\\left(j,n\\right)-O\\left(j,n\\right)\\). The diffusivity coefficients are calculated as harmonic means of values from the adjacent cells. Equation [(2.25.15)](#equation-24-15) is solved for gaseous and aqueous concentrations above and below the water table, respectively. The _R_ term ensure the total mass balance in both phases is properly accounted for. An analogous relationship can be generated for the fully implicit case by replacing _n_ by _n+1_ on the _C_ and _S_ terms of equation [(2.25.15)](#equation-24-15). Using an average of the fully implicit and fully explicit relationships gives:

(2.25.16)[¶](#equation-24-16 "Permalink to this equation")\\\[\\begin{split}\\begin{array}{l} {-\\frac{1}{2\\Delta x\_{j} } \\frac{D\_{m1}^{} }{\\Delta x\_{m1}^{} } C\_{j-1}^{n+1} +\\left\[\\frac{R\_{j}^{n+1} }{\\Delta t} +\\frac{1}{2\\Delta x\_{j} } \\left(\\frac{D\_{p1}^{} }{\\Delta x\_{p1}^{} } +\\frac{D\_{m1}^{} }{\\Delta x\_{m1}^{} } \\right)\\right\]C\_{j}^{n+1} -\\frac{1}{2\\Delta x\_{j} } \\frac{D\_{p1}^{} }{\\Delta x\_{p1}^{} } C\_{j+1}^{n+1} =} \\\\ {\\frac{R\_{j}^{n} }{\\Delta t} +\\frac{1}{2\\Delta x\_{j} } \\left\[\\frac{D\_{p1}^{} }{\\Delta x\_{p1}^{} } \\left(C\_{j+1}^{n} -C\_{j}^{n} \\right)-\\frac{D\_{m1}^{} }{\\Delta x\_{m1}^{} } \\left(C\_{j}^{n} -C\_{j-1}^{n} \\right)\\right\]+\\frac{1}{2} \\left\[S\_{j}^{n} +S\_{j}^{n+1} \\right\]} \\end{array},\\end{split}\\\]

Equation [(2.25.16)](#equation-24-16) is solved with a standard tridiagonal solver, i.e.:

(2.25.17)[¶](#equation-24-17 "Permalink to this equation")\\\[aC\_{j-1}^{n+1} +bC\_{j}^{n+1} +cC\_{j+1}^{n+1} =r,\\\]

with coefficients specified in equation [(2.25.16)](#equation-24-16).

Two methane balance checks are performed at each timestep to insure that the diffusion solution and the time-varying aggregation over inundated and non-inundated areas strictly conserves methane molecules (except for production minus consumption) and carbon atoms.

