Hello,
I had written an answer but it seems lost.
In fact this has been done in Telemac-3D but not in Telemac-2D, though it is certainly possible. What is a boundary term in Telemac-3D would become a source-sink term at all points in 2D. In your case you would need two tracers, namely your specific tracer, and the water temperature.
A user has recently developed the effect of evaporation on tracers, which exists also in Telemac-3D.
The only problem is that in 2D any exchange at the free surface is immediately averaged on the depth, this is a strong and debatable assumption, for example thermal stratifications will not be seen, and the results may be so so...
With best regards,
Jean-Michel Hervouet