## 2 Theoretical Model for Water Temperature Field
### Gas-Liquid Two-Phase Interface Heat Flux Calculation Model
The standard volumetric expressions for the two fluxes and stress components are as follows:  
\[\psi_{s}=\rho_{a}c_{pa}C_{h}U(T_{sw}-T_{sa}) \tag{5}\]  
\[\psi_{l}=\rho_{a}L_{c}C_{e}U(q_{sa}-q_{sa}) \tag{6}\]  
\[\tau=\rho_{a}C_{d}U(u_{si}-u_{i}) \tag{7}\]  
where \(C_{h}\) represents the sensible heat transfer coefficient, \(U\) denotes the mean wind speed at height \(z_{r}\) relative to the water surface, \(T_{sw}\) is the water surface temperature, and \(T_{sa}\) is the air temperature at height \(z_{r}\); \(C_{e}\) signifies the latent heat transfer coefficient, \(q_{sa}\) is the specific humidity at the air-water interface, and \(q_{ta}\) is the specific humidity at height \(z_{r}\); \(C_{d}\) is the stress transfer coefficient, \(u_{si}\) is the water surface velocity, and \(u_{i}\) is the horizontal wind speed component relative to a fixed coordinate system. To determine these transfer coefficients, they are decomposed into relatively independent neutral transfer coefficient components, which are expressed as follows:  
\[ \begin{array}{l}C_{h}=c_{t}^{1/2}c_{d}^{1/2}\\ C_{e}=c_{q}^{1/2}c_{d}^{1/2}\\ C_{d}=c_{d}^{1/2}c_{d}^{1/2}\end{array}  \tag{8}\]  
where the neutral transfer coefficients \(c_{T}^{1/2}\), \(c_{q}^{1/2}\), and \(c_{d}^{1/2}\) are reformulated using the flux functions derived from MOS theory, which are expressed as follows:  
\[ \begin{array}{l}c_{T}^{1/2}=c_{Ttn}^{1/2}/[1-\frac{c_{Tn}^{1/2}}{a_{n}c_{h} }(5)]\\ c_{q}^{1/2}=c_{qtn}^{1/2}/[1-\frac{c_{qn}^{1/2}}{a_{n}c_{h}}(5)]\\ c_{d}^{1/2}=c_{d}^{1/2}/[1-\frac{c_{qn}^{1/2}}{a_{n}c_{h}}(5)]\end{array}  \tag{9}\]where \(f_{k}\) is the von Karman constant, taken as 0.4 [2]; \(a_{n}\) is the scalar correction factor, taken as 1.0 [38]; and \(\zeta\) is the MOS stability parameter, satisfying the following equation:  
\[\xi=z_{r}/L \tag{10}\]  
\[L^{-1}=\frac{f_{k}g}{T_{sa}}(T_{*}+0.61T_{sa}q_{*})/u_{*}^{2} \tag{11}\]  
where \(L\) is the Obukhov length, g is the gravitational acceleration, and \(\phi\) is the MOS profile function incorporating the standard Kansas-type model, satisfying the following equation [39]:  
\[\phi_{h}=\frac{1}{1+\xi^{2}}\phi_{hk}+\frac{\xi^{2}}{1+\xi^{2}}\phi_{c} \tag{12}\]  
\[\phi_{u}=\frac{1}{1+\xi^{2}}\phi_{uk}+\frac{\xi^{2}}{1+\xi^{2}}\phi_{c} \tag{13}\]  
where  
\[\phi_{hk}=2\ln(\frac{1+(1+\gamma_{1}\xi)^{1/2}}{2}) \tag{14}\]  
\[\phi_{uk}=2\ln(\frac{1+(1-\gamma_{2}\xi)^{1/4}}{2})+2\ln(\frac{1+(1-\gamma_{2} \xi)^{1/2}}{2})-2\arctan((1-\gamma_{2}\xi)^{1/4})+\frac{\pi}{2} \tag{15}\]  
\[\phi_{c}=1.5\ln(\frac{(1-\gamma_{3}\xi)^{2/3}+(1-\gamma_{3}\xi)+1}{3})-\sqrt{ 3}\text{arctan}(\frac{2(1-\gamma_{3}\xi)+1}{\sqrt{3}})+\frac{\pi}{\sqrt{3}} \tag{16}\]  
where \(\gamma_{1}\), \(\gamma_{2}\), and \(\gamma_{3}\) are empirical constants. In Equation (9), the subscript \(n\) represents the specific scenario when \(\phi=0\). In this case, the neutral transfer coefficients are directly associated with their corresponding roughness lengths, which are expressed as follows:  
\[ \begin{array}{l}c_{Tn}^{1/2}=\frac{a_{n}f_{k}}{\log(z_{r}/z_{oT})}\\ c_{qn}^{1/2}=\frac{a_{n}f_{k}}{\log(z_{r}/z_{oT})}\\ c_{dn}^{1/2}=\frac{f_{k}}{\log(z_{r}/z_{o})}\end{array}  \tag{17}\]  
where \(z_{o}\), \(z_{oq}\), and \(z_{oT}\) represent the roughness lengths for velocity, humidity, and temperature, respectively. These roughness lengths quantify the strength of heat and moisture transfer via molecular diffusion processes within the thin layer near the water surface. The velocity roughness length satisfies the following equation:  
\[z_{o}=0.011\frac{u_{*}^{2}}{g}+0.11\frac{v_{a}}{u_{*}} \tag{18}\]  
where \(v_{a}\) represents the kinematic viscosity of air.