The time integration of 1D thermal conduction equation of grand surface (equation (55) in Part I is performed by the Crank-Nicolson scheme. The space differencing is evaluated by the second order centered scheme. The grand temperature and depth are evaluated on the grid point and the heat flux is evaluated on the half grid point. The number of vertical grid point is and the suffix of the lowest grid point is . The is assumed to the surface temperature . The finite difference 1D thermal conduction equation is represented as follows.

(66) |

or,

(67) |

where . This equation can be represented in matrix form as follows.

where . The elements of , are represented as follows.

Considering the boundary condition of upper and lower boundaries,
(68) is modified as follows.

(69) |

(70) |

where the elements of and are modified as follows.

is a column vector whose dimension is
are represented as follows.