Can surface sensible heat flux be calculated using the equations from the IFS documentation?


EDIT: The error of the surface sensible heat flux being > 2 orders of magnitude too large is solved and was due to me foolishly neglecting brackets around the denominator in the equation for CH.

It thus seems like the flux can at least be somewhat estimated this way (although the correlation with the flux from the output is still surprisingly low, especially for cases with stable stratification).


Dear all,

I am trying to calculate the surface sensible heat flux from ERA5. However, although my resulting flux has relatively similar spatial patterns, it is more than two orders of magnitude too large, compared to the heat flux directly available in the output (the reason I’m not using the fluxes from the output is because in a later step, I want to calculate fluxes for fixed surface temperatures, to separate atmospheric from surface influences). Of course the error(s) can result from many different things, but when I change the variables slightly to see what happens, I always end up with a largely similar (ca 2 orders of magnitude too large) flux. This seems to suggest that it might be an error related to units, or, that my method is rubbish and can’t be applied in this way. 

I have data for one time step, and I’m checking my flux against the instantaneous surface sensible heat flux (ishf), which is given in W/m^2. As far as I understand, since I calculate my flux from temperature, wind, etc at that particular time step, J and W should be the same anyway, since W =J/s, and the time step is 1 s. Or do I somehow need to scale my resulting flux?

Or, is it simply not possible to calculate the flux using these methods outside of IFS, using the ERA5 data as input?

I would be very happy if someone had any idea about this!

Further details are found below.

Best regards,


I am using the equations from the IFS Documentation CY48R1 - Part IV: Physical Processes

to derive the surface sensible heat flux H (Equation 8.6) from wind speed, the temperature difference between skin temperature and temperature at 2 m, some constants, and the heat transfer coefficient CH. I derive CH using Equation 3.17, with the results from the stability functions psiM and psiH for unstable (Equation 3.22) and stable (Equation 3.24) conditions, respectively.

To get the stability parameter zeta, I have calculated the inverse Obukhov length 1/L following the instructions here:

And since zeta is defined as z/L, 

zeta = z*1/L

where z = 2, since the temperature is valid at 2 m height


For unstable conditions (zeta <0): Equation 3.22

x = (1-16*zeta)^1/4

psiM = pi/2 -2*atan(x)+log((1+x)^2*(1+x^2)/8)

psiH = 2*log(1+x^2/2)

where log is the natural logarithm ln

For stable conditions (zeta >0): Equation 3.24

a = 1, b = 2/3, c=5, d = 0.35

psiM = -b*(zeta-(c/d))*exp(-d*zeta)-a*zeta-b*c/d

psiH = -b*(zeta-(c/d))*exp(-d*zeta)-(1+(2/3)*a*zeta)^1.5-b*c/d+1


k = 0.4; 

k^2 /( log((z+z0M)/z0M)-psiM*((z+z0M)/L)+psiM(z0M/L) )* ...

( log((z+z0M)/z0H)-psiH((z+z0M)/L)+psiH(z0H/L) )

where z0M is the forecast surface roughness (fsr)

and z0H is the surface roughness for heat, which I derived from "Forecast logarithm of surface roughness for heat (flsr)" as z0H= e^flsr


H = rho * cp * U * CH * (T +gz/cp - T_sk)


g = 9.80665

cp = 1004.709

rho is air density, approximately 1.29

U is wind speed

T is temperature at 2 m

T_sk is the skin temperature