If the temperature of the pipe wall (i.e., the surrounding rock) can be considered a constant, then the existing EPANET program and interface can be used by making the following change to the wallrate() function at the bottom of the quality.c module: replace the code segment underneath the "if (WallOrder == 0.0)" with the following:

if (WallOrder == 0.0)

{

return kf * 4.0/d * (Climit - c)

}

Then re-compile the epanet2.dll library. When you run EPANET, set the
Wall Reaction Order to 0, Global Bulk reaction coefficient to 0, the
Global Wall coefficient to any non-zero value, the Limiting
Concentration to the constant pipe wall temperature, and the Relative
Diffusivity to the thermal diffusivity of water (in sq ft /day) divided
by 0.00112.

These changes express the rate of temperature change in the bulk water as the product of a heat transfer coefficient (kf, computed internally by EPANET), the area per unit volume of the pipe (4/d), and the difference between the wall and water temperatures (Climit - C). If the heat transfer mechanism is more complicated than this then it might be better to follow Steve's suggestion to use the MSX toolkit.

We noticed that in epanet22, the quality.c had splitted into 3 files.
Now we need to make the above adjustment in qualreact.c, line 431:

    if (qual->WallOrder == 0.0)             // 0-order reaction */

    {

        kf = SGN(kw) * c * kf;              //* Mass transfer rate (mass/ft2/sec)

        kw = kw * SQR(pr->Ucf[ELEV]);       // Reaction rate (mass/ft2/sec)

        if (fabs(kf) < fabs(kw)) kw = kf;   // Reaction mass transfer limited

        return (kw * 4.0 / d * (qual->Climit -c));              // Reaction rate (mass/ft3/sec)

    }