How to diagnose TSK and T2 in an urban grid cell?

TSK and T2 are the most widely examined variables, yet their calculations are quite convoluted. Let's start with the basics.

Let's refrain ourselves in the world of WRF-Noah (without the mosaic approach)-SLUCM (single-layer urban canopy model). In this case, a grid cell, if the dominant land cover is urban, will be separated into an impervious part (with the fraction of FRC_URB2D) and a pervious part (with the fraction of 1-FRC_URB2D). The impervious part is handled by SLUCM, and the pervious part is handled by Noah.

The first thing to note is that TSK may be a prognostic variable while T2 is almost always a diagnostic variable. A prognostic variable means that it has to be computed at every single time step, otherwise the model crashes. On the other hand, a diagnostic variable is something nice to have, but not necessary.

In Noah, TSK is indeed a prognostic variable so it will be always computed for the pervious part. Let's call it T1 as in module_sf_noahdrv.F.

Its counterpart TSK_urban, from the SLUCM, is however a diagnostic variable. Wait, how is this possible, shouldn't the model need TSK_urban at each time step to move forward? Well, not really. The SLUCM further separates the impervious part of urban land into roof, wall, and ground (let's call them 3 facets) and solves the surface energy balance for each facet, yielding a TSK for each facet. The TSK_urban is then diagnosed based on the surface energy budgets and states of all facets.

One could come up with different ways of diagnosing TSK_urban but the SLUCM diagnoses TSK_urban based on the bulk sensible heat flux formulation, with the basic premise being that using TSK_urban should give the same total sensible heat flux from the impervious urban land. While there is nothing wrong with this approach, one has to be careful with the aerodynamic resistance ($r_a$) or the equivalent conductance used in the bulk sensible heat flux formulation, see below. It's clear that with everything else being the same, a larger $r_a$ leads to a larger $T_s$, in this case, a larger TSK_urban.

$$ H = \rho c_p \frac{T_s - \theta_a}{r_a} $$

Now comes the most intriguing part of the SLUCM code: the $r_a$ is calculated outside the SLUCM code loop using the roughness lengths of the pervious land (note that since $r_a$ is a resistance for heat transfer instead of momentum transfer, it will depend on both the momentum and thermal roughness lengths). This is rather bizarre because $r_a$ should be, at least, related to the impervious land property like building height, canyon width, etc. Yet, in SLUCM $r_a$ has nothing to do with the impervious part of urban land. Why is that?

Note that in the model we can use whatever $r_a$ value for diagnosing TSK_urban, because again TSK_urban is a diagnostic variable that does not influence any following calculations. This should be contrasted with the prognostic equations for TSK of roof, wall, ground, in which the $r_a$ was using the properties of the impervious part of urban land. So using the $r_a$ computed with The roughness lengths of the pervious land to diagnose TSK_urban is not wrong per se, but this decision will have implications for how we interpret TSK_urban.

Now that we understand that TSK_urban is a diagnostic variable and how it is diagnosed, let's see how is the TSK of the entire grid cell is diagnosed. It is quite simple actually. The TSK of the entire grid cell is diagnosed as

TSK = FRC_URB TSK_urban + (1-FRC_URB) T1

Let's use new notations, with subscripts u and r indicating impervious and pervious respectively. $f_u$ is the impervious fraction. So the TSK of the entire grid cell is:

$$ T_s = f_u T_s^u + (1 - f_u) T_s^r $$

The total sensible heat flux of the grid cell is formulated in a similar way:

$$ H = f_u H^u + (1 - f_u) H^r = f_u \rho c_p \frac{T_s^u - \theta_a}{r_a^u} + (1 - f_u) \rho c_p \frac{T_s^r - \theta_a}{r_a^r} $$

Now comes the punch line: the above equation can be simplified if $r_a^u = r_a ^r = r_a$

$$ H = f_u H^u + (1 - f_u) H^r = f_u \rho c_p \frac{T_s^u - \theta_a}{r_a^u} + (1 - f_u) \rho c_p \frac{T_s^r - \theta_a}{r_a^r} = \rho c_p \frac{f_u T_s^u + (1 - f_u) T_s^r - \theta_a}{r_a} = \rho c_p \frac{T_s - \theta_a}{r_a} $$

The above equation is nice in the sense that the grid cell sensible heat flux $H$ can be related to the grid cell TSK ($T_s$) through $\rho c_p \frac{T_s - \theta_a}{r_a}$. This is particularly important for setting the stage to diagnose T2, as you will see later. But this comes with the price that we have to set $r_a^u = r_a^r$. This is perhaps why in the SLUCM code, the calculation of $r_a^u$ is all commented out and the code ends up using $r_a^r$, which is computed outside the SLUCM code and computing using roughness length of the pervious part of urban land.

Now the question is, yes, maybe we should set $r_a^u = r_a^r$. But in this line of argument, $r_a^u$ and $r_a^r$ values can be anything as long as they are equal. Why do we need to set their values to be those computed using roughness length of the pervious part of urban land?

This is certainly a valid question and in my opinion, using the $r_a$ computing only using roughness length of the pervious part of urban land is wrong. Again, $r_a$ should, at least, reflect the properties of the impervious part of urban land. For example, one way to address this issue is to use a composite roughness length based on the roughness lengths of both pervious and impervious land.

To do so, you need to change ZNT from a IN variable to be a INOUT variable for the module_sf_urban.F, which I did 10 years ago when I implemented the mosaic approach. Now we have roughness length of both pervious and impervious land, one way to composite them is to add the following line of code in module_sf_noahdrv.F after CALL urban:

        ZNT(I,J)= EXP(FRC_URB2D(I,J)*ALOG(ZNT_URB)+(1-FRC_URB2D(I,J))* ALOG(Z0K))

This change has effects that are reasonably predictable. As the roughness length of the pervious land is usually smaller than the roughness length of the impervious land, the composited roughness length is then often larger than the roughness length of the pervious land. With everything else being equal, the TSK will be reduced.

Below are two simulations for Boston. Just looking at the left panels which show the TSK of the entire grid cell. The 1st one is from the original code, and the 2nd one is from a modified code with the proposed change. You can see that the TSK drops. I also did two other simulations where the roof albedo is increased from 0.2 to 0.9 (mid-panels) to quantify the effects of white roofs (right panels). You can also see that the effects of white roofs are reduced. Again, not too surprising because the proposed code change essentially introduced more mixing.

TSK_test1

TSK_test2

Now that we understand the effects of the proposed code changes. What does it represent physically? Well, the two scenarios represent two extremes. In the original code, the mixing over the pervious land is only produced by the roughness elements of the pervious land and is completely independent of the imperious land. You can think of this as a scenario where the pervious part is a park far away from buildings. In the proposed change, the mixing over the pervious land is also affected by roughness elements over the impervious land. You can think of this as a scenario where the pervious part is grass dispersed around buildings.

I would argue that only the second (or the proposed) approach is internally consistent with the assumption we have made earlier on, that is $r_a^u = r_a^r$. What $r_a^u = r_a^r$ implies is that the pervious and impervious land feels the same mixing power of turbulence, meaning that the pervious and impervious land are so intertwined/dispersed within each other. Then then only consistent way of computing this mixing is to use the roughness lengths of both pervious and impervious land.

Good. Now are you ready for T2? In WRF, T2 is diagnosed as

$$ H = \rho c_p \frac{T_s - \theta_a}{r_a} = \rho c_p \frac{T_s - T2}{r{a,2}} $$

Can you now see why we need to assume $r_a^u = r_a^r = r_a$ when we diagnose the grid cell $T_s$ or TSK? Because only by doing so can we write the grid cell $H$ as $\rho c_p \frac{T_s - \theta_a}{r_a} = \rho c_p \frac{T_s - T2}{r{a,2}}$ where $r_{a,2}$ is the aerodynamic resistance between the surface and 2 meters about the surface (we haven't defined surface yet, have we?). Using the above equation we can obtain the grid cell $T_2$ as

$$ T_2 = Ts - \frac{H r{a,2}} {\rho c_p} $$

Incidentally, in the original code $r{a,2}$ is also computed using the roughness length of the pervious land. With the proposed modification, $r{a,2}$ is computed using a composite roughness length. Let's see how the results have changed. Similarly, the T2 values are reduced and the white roof effects are reduced. This again can be explained by the stronger mixing introduced by the modified code.

T2_test1

T2_test2

In essence, WRF-Noah-SLUCM tries to diagnose a grid cell TSK first, and then a grid cell T2. Both are diagnosed based on the premise that the grid cell sensible heat flux can be computed using the grid cell TSK or the grid cell T2. But (I think) the calculation of $ra$ and $r{a,2}$ in WRF-Noah-SLUCM is not internally consistent with the diagnosis of grid cell TSK and T2 unless we adopt the modification proposed here.

A broader and remaining question is: is there a better way to diagnose TSK and T2?

DanLi-BU / WRF

How to diagnose TSK and T2 in an urban grid cell? #5