New option for SLUCM to consider global distributed urban parameters

[epn09]

New option for SLUCM to consider global distributed urban parameters

TYPE: new feature

KEYWORDS: SLUCM, urban parameters, anthropogenic heat

SOURCE: Do Ngoc Khanh (Tokyo Institute of Technology)

DESCRIPTION OF CHANGES:
This PR adds a new feature to WRF SLUCM by allowing consideration of spatially varying global distributed urban parameters and spatially hourly monthly varying anthropogenic heat.

ISSUE: N/A

LIST OF MODIFIED FILES:
M Registry/Registry.EM_COMMON
M Registry/registry.dimspec
M dyn_em/module_first_rk_step_part1.F
M dyn_em/module_initialize_real.F
M phys/module_pbl_driver.F
M phys/module_physics_init.F
M phys/module_sf_clm.F
M phys/module_sf_noahdrv.F
M phys/module_sf_urban.F
M phys/module_surface_driver.F
M phys/noahmp
M share/output_wrf.F

TESTS CONDUCTED:

1. The modification has been tested and used in previous publications.

• Initial development: Varquez, A. C. G., Nakayoshi, M., & Kanda, M. (2015). The effects of highly detailed urban roughness parameters on a sea-breeze numerical simulation. Boundary-layer meteorology, 154, 449-469. https://doi.org/10.1007/s10546-014-9985-4
• Global extension: Khanh, D. N., Varquez, A. C., & Kanda, M. (2023). Impact of urbanization on exposure to extreme warming in megacities. Heliyon, 9, e15511. https://doi.org/10.1016/j.heliyon.2023.e15511

3. Are the Jenkins tests all passing?

RELEASE NOTE: This modification add a two options (use_distributed_aerodynamics and distributed_ahe_opt) to WRF SLUCM (sf_urban_physics = 1) so that spatially varying urban morphological parameters (building height, plan area index, frontal area index, roughness length for momentum, and displacement height) can be considered.

[epn09]

This PR goes in pair with NCAR/noahmp#96. Without it, the code will not compile.

[weiwangncar]

@epn09 The compilation for the code has failed.
@cenlinhe Where is 'noahmp#96'?

[cenlinhe]

I think @epn09 also modifies the Noah-MP driver part of the code, but the way of submitting this PR to include the Noah-MP submodule update may not be correct, which may be the cause for the regression test to fail.

[cenlinhe]

I will come back to this PR later for a more detailed review when I have more time.

[cenlinhe]

@tslin2, you may also want to review this urban code which is also related to Noah-MP PR I assigned you as a reviewer. Note that this is just for Noah-MP v4.5 not the latest refactored code.

[tslin2]

@tslin2, you may also want to review this urban code which is also related to Noah-MP PR I assigned you as a reviewer. Note that this is just for Noah-MP v4.5 not the latest refactored code.

okay, will take a look.

[epn09]

@cenlinhe May I ask if there's any update?

[cenlinhe]

This PR requires an update from Noah-MP driver part for v4.5 first and then connect the WRF update with the updated Noah-MP commit from Noah-MP GitHub in a new PR. I will need to update the Noah-MP GitHub commit first and then probably have to re-submit a PR with the updated Noah-MP commit for this.

[cenlinhe]

@epn09 what is the difference between LF_URB2D_S and LF_URB2D? Why do you need to create this new variable instead of using the existing LF_URB2D?

[epn09]

@cenlinhe S in LF_URB2D_S means "single". If my understanding is right, LF_URB2D has 4 values for 4 wind direction, but our dataset only have one value for all wind direction. So, I created this new variable.

[cenlinhe]

OK, I see. This makes sense.

[cenlinhe]

@dudhia What is the difference between "i" and "i01"?

[weiwangncar]

@cenlinhe Here letter i represents input data. The number 0 means it is part of wrfinput, output from real, and input to the model, and 1 is the input to real program or the met_em output data.

[cenlinhe]

@weiwangncar If we do not specify "0" and/or "1" after "i", does it mean it can be in any input files or by default assume it is from wrfinput? For example, for variables like HGT_URB2D, they are only specified as "i" instead of "i01".

[cenlinhe]

@weiwangncar @dudhia Since these variables are optional (i.e., users may not provide these in their geo_em input file), I think it should be set to "i" instead of "i01". Otherwise, WRF will look for these variables in wrfinput or met_em files, and report errors? Is this correct?

[weiwangncar]

@cenlinhe Sorry I missed these comments. If you just have an 'i', that variable is only appearing in wrfinput - whether it has value will depend on what real does. If there is no input to real, the array may be all zero. If the data comes from geogrid and metgrid, it should be set to i01.

[cenlinhe]

need to add config_flags%sf_urban_physics == 1 in the if-statement too.

[cenlinhe]

need to also include a treatment for urban LCZ type so that the model will not produce weird values for LCZ pixels. If this scheme does not work with LCZ, then add a if-statement here to show error message for this incompatibility.

[epn09]

@cenlinhe This is fixed.

[epn09]

@cenlinhe I added if-statement and error messages to share/module_check_a_mundo.F.

[cenlinhe]

@epn09 not sure if this is a good way to estimate urban fraction (=1-vegfrac), why not using the LAND_USE_F variable to read in the urban fraction? Also, I think the shdmax should be used instead of shdavg.

[epn09]

@cenlinhe Please correct me if I'm wrong but LANDUSEF is not so good because it's derived from land use map (for example, MODIS). If a WRF grid will all MODIS grids classified as urban, the WRF grid will have LANDUSEF(urban) = 100. But from our experience and from reality, 100% urban cover is really rare. So we assume that urban fraction + vegetation fraction = 1, because even in 100% urban cover grid there is still some vegetation fraction in the input data.

As for the shdmax, I think if we do that urban fraction might be a little bit small because trees in its growth season can cover very high fraction of the ground. On the other hand, shdmin may make urban fraction too big. Hence, I use the value in the middle shdavg.

[cenlinhe]

OK, this makes sense, at least for now. more tests with the impact of shdavg vs shdmax are needed in the future.

[DanLi-BU]

@tslin2 @epn09 OK, let me comment here since this is probably the first place where it showed up. I think I understood what Khanh is trying to do here. If there is no other data source for FRC_URB2D, one can come up with a few ways to assign FRC_URB2D values. It used to be that FRC_URB2D is prescribed in the URBPARM.TBL (e.g., 0.5 for low intensity residential urban etc). With NLCD data, FRC_URB2D is actually estimated from NLCD and is no longer using prescribed values from URBPARM.TBL. Now this part of code is trying to assign FRC_URB2D based on vegetation fraction, which is not a terrible idea. But this will limit the compatibility of this code with potentially new land use/land cover data that can provide FRC_URB2D.

Khanh (first of all, congratulations on becoming a doctor, way to go), regarding If a WRF grid will all MODIS grids classified as urban, the WRF grid will have LANDUSEF(urban) = 100, this is not true actually. LANDUSEF is always computed based on land use land cover data in WPS, and stored and passed to WRF. But it is not used when the mosaic approach is not enacted. In a somewhat similar way, FRC_URB2D can be also computed using land use land cover data (especially with high resolution data like NLCD) and it is actually used by WRF (at least before LCZ was introduced, I don't know the details after LCZ is introduced).

[epn09]

@DanLi-BU Thanks for your note about NLCD. Well-noted.

[cenlinhe]

@epn09 This may not be the best place to add this anthropogenic heating flux. First, if putting it in this pbl driver module, should it be added to a place after land surface model and urban scheme is called? Second, it may be better to include this to the SLUCM code module instead of pbl driver to avoid messing up with other physics.

[epn09]

@cenlinhe The reason we didn't add AH flux to the urban scheme is because the construction of the AH dataset did not assume that AH comes from only urban.

• Energy loss are distributed equally for all grids within a country.
• There are strong point sources of AH (power station, steel mills) that do not necessarily belong to urban grids.

What's your opinion on this?

The full paper is available here (https://www.nature.com/articles/s41597-021-00850-w).

[cenlinhe]

OK, I see. This makes sense to me, as the AH is not only for urban AH. power station may also belong to urban?

[epn09]

@cenlinhe Power station may be in either rural or urban. When constructing the AH dataset we did not make that distinction.

[cenlinhe]

double check consistency for the vegfrac used in urban fraction (=1-vegfrac) and the original vegfra used in Noah.

[cenlinhe]

should this already be done above in the initialization and real part?

[epn09]

@cenlinhe To be honest I also don't know why I must do this. But as you can see, the code just below this line for wudapt also do the same thing for FRC_URB2D so I just follow the way the existing code was written.

                  IF(UTYPE_URB==1) FRC_URB2D(I,J)=1.
                  IF(UTYPE_URB==2) FRC_URB2D(I,J)=0.99
                  IF(UTYPE_URB==3) FRC_URB2D(I,J)=1.00
                  IF(UTYPE_URB==4) FRC_URB2D(I,J)=0.65
                  IF(UTYPE_URB==5) FRC_URB2D(I,J)=0.7
                  IF(UTYPE_URB==6) FRC_URB2D(I,J)=0.65
                  IF(UTYPE_URB==7) FRC_URB2D(I,J)=0.3
                  IF(UTYPE_URB==8) FRC_URB2D(I,J)=0.85
                  IF(UTYPE_URB==9) FRC_URB2D(I,J)=0.3
                  IF(UTYPE_URB==10) FRC_URB2D(I,J)=0.55
                  IF(UTYPE_URB==11) FRC_URB2D(I,J)=1.

[tslin2]

this part is the issue in mosaic where the 2d spatially urban fraction is not used. So after the first timestep, the urban fraction is different compared to the initial 2d data. It should already be done in the urban_init_var

[cenlinhe]

I will need to take a more detailed look at this. These FRC_URB2D assignments for each LCZ are supposed to be reading from the URBPARM_LCZ.TBL table values instead of hard-coded.

[cenlinhe]

@epn09 @tslin2 This original part of re-assigning FRC_URB2D is wrong. FRC_URB2D is initialized in module_physics_init.F and should only be read in as a intent(IN) variable without modifying. This has been done in both Noah and NoahMP driver. So this part should be removed from lsm_mosaic module.
@tslin2 you may want to test if removing this would solve your mosaic error.

[tslin2]

maybe @DanLi-BU can also comment on this.

Thanks

[DanLi-BU]

I haven't read the WRF code this closely for such a long time but it's good to be back in the coding world. I don't think I added this line of FRC_URB2D(I, J) = MAX(0.1, MIN(0.9, 1 - shdavg(I, J) / 100.)) when I implemented the mosaic and I agree FRC_URB2D should not be changed (unless there is a strong justification).

[tslin2]

FRC_URB2D(I, J) = MAX(0.1, MIN(0.9, 1 - shdavg(I, J) / 100.))

is proposed by @epn09 for new option to consider global distributed urban parameters. @cenlinhe and I are reviewing.

[tslin2]

@DanLi-BU you can also review if you like

[DanLi-BU]

Khanh has to buy me a lot of coffee for me to review the entire code :) I'll try my best to comment here and there but I can't really promise.

[cenlinhe]

also need to add a check for cases where distributed_aerodynamics_option is true but the input variables/data are not provided (i.e., MH_URB <=0 or LP_URB<=0).

[epn09]

@cenlinhe Within this if statement, some default values for the parameters when they are not available are set.

ZR = MAX(MH_URB, 3.5)
Z0C = MAX(Z0_URB, 0.1)
lambda_p = MAX(0.05, MIN(1.0, LP_URB))
lambda_f = MAX(0.05, MIN(1.0, LF_URB_S))

[cenlinhe]

OK, I see.

[cenlinhe]

@epn09 It seems you replaced original TCP related quantity with TA, why is this? based on assumptions for your distributed parameter generation?

[epn09]

@cenlinhe According to the materials that I read, in this scheme, canopy temperature is not calculated (hence, TCP and TC are unused). Heat from all surfaces (roof, wall, roads) are released to the upper atmosphere directly. Therefore, TA is used.

References:

• Section I of http://www.ide.titech.ac.jp/~kandalab/download/WRF_URBAN/add_improvements.html
• https://link.springer.com/article/10.1007/s10546-004-7956-x

[cenlinhe]

@epn09 not sure if this is a good assumption. maybe still use the logscale decay? Is this assumption related to how you treat the roughness from your new input parameters?

[epn09]

@cenlinhe Currently in this parameter database that we have, there is no grid with z0 higher than 10 m, so this if statement actually will never be executed. But the treatment of grid with z0 > 10 m (if exists) is something we also don't know how to properly handle because in my understanding log-law is not valid for z < z0. So if z = 10 m and z0 > 10 m, I don't really know what should be done. In the mean time, this if statement can be removed. Please let me know your opinion.

[cenlinhe]

I do not have a good solution to this either. The default treatment is using logscale decay. For Z0C > 10, the log law should give a relatively small value (close to 0) at z=10. Maybe for now, just follow the default treatment.

[epn09]

@cenlinhe This is fixed. The if Z0C > 10 is removed.

[cenlinhe]

@epn09 your distributed parameters are only for urban type = 2?

[epn09]

@cenlinhe In my understanding, the difference between urban types (low intensity, high intensity, LCZ 1 to 10) in WRF right now are only the differences in the parameters in URBPARM.TBL. When we input a spatial map of these urban parameters from outside, the value of urban type doesn't matter any more. So, it's just set to 2 for no special reasons.

[cenlinhe]

OK

[cenlinhe]

I have reviewed this PR and have provided a few comments and questions in some code lines. This new capability of reading in and using global spatially varying distributed urban parameters are very important. However, it seems that the current implementation in this PR has a few limitations (correct me if I am wrong):
(1) requires specific input urban parameter data from the developer's group. Did the develop also submit a PR to WPS to include their data in the WRF system? Otherwise, it is very hard for other users to use this capability.
(2) The current implementation seems incompatible with the existing I/O capability from using NUDAPT parameter data. However, these two I/O systems should be similar but of course reading in different datasets and variables. Is there a way to synthesize these two to be more general?
(3) This update does not work with LCZ urban type, which may be OK for now.
(4) This update seems only be applicable for traditional urban type = 2 based on the comment I provided in one line above (maybe I misunderstood it). This strongly limits the code application.

[epn09]

@cenlinhe Thank you very much for reviewing. Some of your comments need code modification so I'll address them separately. But below is my response for the 4 points you raised.

(1) requires specific input urban parameter data from the developer's group. Did the develop also submit a PR to WPS to include their data in the WRF system? Otherwise, it is very hard for other users to use this capability.

No, I haven't submit a PR. Please let me know the procedure and I'll do it promptly.

(2) The current implementation seems incompatible with the existing I/O capability from using NUDAPT parameter data. However, these two I/O systems should be similar but of course reading in different datasets and variables. Is there a way to synthesize these two to be more general?

In my understanding, NUDAPT only change the urban parameters but do not change the physics. However, our scheme not only change the urban parameters but also represent the urban canopy as a whole and disregard the individual effects of the roof and canyon. More details can be found here (http://www.ide.titech.ac.jp/~kandalab/download/WRF_URBAN/add_improvements.html). Please ignore section II of the link because it's obsolete thanks to the implementation of WRF mosaic option.

So, at the moment, I don't have an idea to synthesize the two approaches.

(3) This update does not work with LCZ urban type, which may be OK for now.

Yes. And I can't think of a valid use case of combining our approach with LCZ for now, even though both aims at the same goal of making more realistic urban representation in WRF.

• LCZ aims at providing more urban categories (from the default 3),
• Our approach uses spatial maps of urban parameters (in other words, there is no urban category anymore).

Of course, the two approach can be compared against each other.

(4) This update seems only be applicable for traditional urban type = 2 based on the comment I provided in one line above (maybe I misunderstood it). This strongly limits the code application.

As mentioned in (3) and replied to your comment. Urban types play no role in this scheme. Type 2 is set in the code for no special reason.

[epn09]

@cenlinhe I think I've addressed all of your comments. Please let me know if you have any further concern.

[cenlinhe]

@epn09 I am going to re-submit a PR with your modifications on urban part and refer back to this PR, because I need to update Noah-MP GitHub with your mods first and re-connect to the correct Noah-MP branch and commit.
In the near future, please submit a PR to WPS GitHub with your updated code that works with your new urban parameter dataset. Please also prepare the WPS binary data format for your dataset and test it in WPS, so that the WRF team could directly incorporate your dataset into WRF/WPS database for users to download.

[cenlinhe]

I have re-submitted this PR at: #1986, which includes the latest develop WRF branch and connect to the new Noah-MP v4.6 release branch.

[weiwangncar]

This has been replaced by PR-1986.