Regarding Dataset prepartion for Kelaniya River Basin, Sri Lanka

[Hrekhani97]

Hello,

I am Madhawa, a student from Sri Lanka. I am planning to use the LISFLOOD model for my research on the Kelaniya River Basin in Sri Lanka. I have successfully installed the model and run the test cases. Currently, I am working on creating a custom dataset for the Kelaniya River Basin. I’ve already developed general maps and topography maps, except for the elvstd.map, and have created land use fraction maps based on local land use data.

I am facing a few challenges and would appreciate your guidance on the following issues:

I used the MERIT DEM to create all the general and topography maps and plan to run the model at a 1 km x 1 km spatial resolution. Could you please advise on how to create the elvstd.map file using PCRaster within QGIS? Do I have to calculate the elevation standard deviation for the original Merit DEM file and resample that to my target resolution?

For the land use fraction maps, I relied on a local land use survey to obtain accurate data. Could you explain the process for calculating land-use-dependent layers, such as the crop coefficient and crop group number? The methodology on the website involves a large number of land use classes, whereas my dataset for the Kelaniya Basin uses a simplified local classification. How should I proceed to develop these two layers in this context?

Thank you very much for your time and assistance. Your support would be greatly appreciated.

Best regards

[StefaniaGrimaldi]

Dear @Hrekhani97,

we are pleased to know that you plan to use OS LISFLOOD in the frame of your research study.

Please find in this page the protocol for the computation of the elvstd map: https://ec-jrc.github.io/lisflood-code/4_Static-Maps_topography/ . The elevation field at the target resolution (1 km) is obtained by computing the average of the original high resolution map (mean.unweighted()). elvstd is the standard deviation of the elevation of all the high resolution cells included into one grid of your target grid (stDev.unweighted()). To achieve this aim, you can use the software and protocol of your choice.

In each pixel, the sum of all land use fractions must be 1. This is a strict requirement to ensure mass balance and correct computations. This page https://ec-jrc.github.io/lisflood-code/4_Static-Maps_land-use-depending/ explains how the classes included in the "Copernicus Global Land Cover Layers: CGLS-LC100 collection 2" were associated to the land cover fractions. When using the local classification available for your study area, you can allocate each class of your database to the most representative land cover fraction (forest, irrigation, other).

As a general piece of advice, we strongly recommend using netcdf format, rather than pcraster format (for example, the netcdf format allows you to define pivotal metadata information such as the geographical reference system).

I hope that this answer was helpful.
Please let me know if I missed or misunderstood parts of your inquiry.
Kind regards,
Stefania

[Hrekhani97]

Dear @StefaniaGrimaldi,

Thank you so much for the information provided. I have now created the land use fraction maps and am currently preparing the land use-dependent maps. I have a quick question regarding the creation of crop coefficient maps such as cropcoef_f and cropcoef_irri. When assigning values, they are specific to the relevant land use classes. For example, the final crop coefficient values for the forest file will only contain pixel values where forest land use classes are present, and other pixels will be assigned no-data values, correct? As you can see, the pixel values are ranging from 0.7-1.1 and no-data value was assigned to the rest of the area (where there is no forest cover). However, I noticed that there are no missing or no-data pixels in your test datasets and the global 03-minute dataset. But, when preparing these files, we assign a 0 value for the areas with no forest cover right, but, how did all the pixels get a value at the end? Could you kindly clarify this?

Below is my forest fraction layer (just for your reference).

Thank you for your assistance

[Hrekhani97]

Hi, I think I figured it out. We need to replace the null values with the global mean values of Kc, Kg, and Km.

Thank you for your assistance.

Best regards

[StefaniaGrimaldi]

Dear @Hrekhani97,
your understanding is correct. OS LISFLOOD does not accept No Data within the modelling domain. As explained here https://ec-jrc.github.io/lisflood-code/4_Static-Maps_land-use-depending/, the recommended solution is to use mean values (computed by excluding all nil Kc, Kg and Km values, respectively) for pixels where the specific land cover class is not present.
We will be pleased to provide further assistance,
kind regards,
on behalf of the developers team,
Stefania

[Hrekhani97]

Dear @StefaniaGrimaldi,

Could you please provide some additional clarification on the following point?

When developing land use-dependent layers, we primarily consider three land use fractions: forest, irrigated crops, and other LULC types. Therefore, when calculating the Manning's coefficient for these three land use classes, we do not take built-up areas/sealed surfaces into account, correct? However, when it comes to flood modeling, the Manning's roughness of sealed surfaces plays a significant role specially during a flood event correct?

If I understand correctly, we are assuming that approximately 25% of the sealed surfaces are permeable. With this in mind, would it be appropriate to categorize this portion under a under other landuse class when defining the manning coefficient?

How do the model accurately incorporate surface roughness details for sealed surfaces in our calculations?

Thank you for your continuous support. I look forward to your guidance on this matter.

[StefaniaGrimaldi]

Dear @Hrekhani97,

OS LISFLOOD computes surface routing for sealed surfaces (included in the direct runoff fraction) in this module: https://github.com/ec-jrc/lisflood-code/blob/fe94c33b6bd9fe09fa7be43a6dd4fdc8a11f744c/src/lisflood/hydrological_modules/surface_routing.py

The Manning's roughness coefficient for sealed surfaces is set equal to 0.02 in this line (this value is applied to all pixels within the computational domain):

lisflood-code/src/lisflood/hydrological_modules/soil.py

Line 140 in fe94c33

self.var.NManning = self.var.defsoil('MapN','MapForestN', 0.02, OrderedDict([self.var.dim_runoff, self.var.dim_pixel]))

The use of a different value or of an input map requires small edits to the source code. In the first case, you will only need to edit the line above by replacing the current value. In the second case (i.e. if you wish to use a map), you will need to prepare a map, add the corresponding key in the .xml settings file, and edit the source code to enable the use of the new input file.

Regarding your question on the permeable fraction of urban areas, you correctly understood the assumption explained in https://ec-jrc.github.io/lisflood-code/4_Static-Maps_land-use/, and your proposed solution is appropriate.

I hope that this answer helps,
kind regards,
on behalf of the developers team,
Stefania

[Hrekhani97]

Dear @StefaniaGrimaldi,
Thank you for your clarification. I am currently preparing the soil hydraulic property layers using the provided pedotransfer functions. I referred to the reference article you mentioned in the methodology (Toth et al. (2015)), and it appears that these PTFs are specific to the European region. In this context, would it be appropriate to use the same PTFs to derive soil hydraulic properties for South Asian countries, such as Sri Lanka?

I have another question: The original bulk density values downloaded from ISRIC are in units of 'cg/cm³'. To apply the PTFs, we need to convert these values to 'g/cm³', correct? Additionally, are there any other soil data layers that require unit conversion before applying the PTFs?

Thank you for your continuous support. I look forward to your guidance on this matter.

[StefaniaGrimaldi]

Dear @Hrekhani97,

thank you for your interesting question. Annex C2 of https://hess.copernicus.org/articles/28/2991/2024/hess-28-2991-2024.html states "Users can decide to derive soil proprieties from different PTFs, but the general principle presented here remains valid". The scientific literature enumerates several pedotransfer functions, as well as comparison studies (here are a just a few examples: https://www.sciencedirect.com/science/article/pii/S0022169419305049, https://www.sciencedirect.com/science/article/pii/S2352009423001165). Modelers can decide to implement pedo-transfer functions which are specific for their area of interest, if available (here is an example of a recent publication focusing on a specific area https://www.sciencedirect.com/science/article/pii/S0016706124002052).

Regarding the units, you are correct, input data must be converted to the units required by the selected pedotransfer functions. This requirement applies to all input data. The conversion of ISRIC bulk density values from cg/cm3 to g/cm3 is correct when applying Toth et. al equations. Another example of input data that requires conversion is the organic carbon content (delivered in dg/kg).

I hope this answer helps,
kind regards,
Stefania

[Hrekhani97]

Dear @StefaniaGrimaldi, @doc78,
Thank you very much for the clarification provided earlier. I have successfully calculated the soil hydraulic parameters for each ISRIC depth [(0-5), (5-15), (15-30), (30-60), (60-100), (100-200)] by applying the PTFs to the corresponding depths of the ISRIC soil characteristics. This step is part of the initial process for calculating the parameters for depth layers 2 and 3. However, I found the next steps a bit confusing, and I would greatly appreciate it if you could clarify a few things for me.

Based on my understanding, the process is as follows:

1. After completing the above step, we need to calculate the weighted average parameter values for soil layers 2 and 3, as the soil depth varies spatially in these layers. In the example provided on the website, the soil depth of the second layer ranges from 5-50 cm. Is this a standard value that should be applied universally, or should this value be adapted according to the specific soil depth of layer 2 in different land use types (e.g., forest, other types)? I have attached an image showing how the soil depths for layers 2 and 3 differ between forest and other land use types in my study area. Specifically, in my case, the soil depth for layer 2 in forest and other LULC types ranges from 0.5 to 2.5 meters, while soil depth for layer 3 ranges from 0.5 to 10.9 meters. Should I still consider 5-50 cm for soil depth layer 2, or should I adjust this to reflect the actual depth ranges in my study area?

2. If I calculate the soil properties as described on the website, how should I combine them with the forest and non-forest fractions? The website suggests calculating the weighted average of the parameter values based on the percentage of forest (or other land use types) per grid cell. However, I am unclear on how this should be done. My understanding is that when multiplying the soil parameters by the forest/other fractions (i.e., the percentage for each cell), some pixels may end up with no values (e.g., 0% forest coverage). How should I address this issue?

I apologize if I am asking too many questions, but I just want to ensure that I am following the correct process.
Thank you for your continuous support. I look forward to your guidance on this matter.

[Hrekhani97]

Dear @StefaniaGrimaldi,
My idea regarding the second question is as follows: Suppose we have a forest fraction layer, which indicates the percentage of forest cover in each grid cell. To calculate the Ksat1_forest layer, we would multiply the Ksat layer by the forest fraction layer. However, there may be some pixels in the forest fraction layer where the value is 0 (indicating no forest cover). In such cases, the corresponding pixels in the Ksat1_forest layer will have a value of 0 for those particular grid cells. Is that correct

[StefaniaGrimaldi]

Dear @Hrekhani97,

please find below our answers:

1. the computation of soil property values must account for the actual depth of the soil layers, according to your maps. You might be interested in the example of computation provided in https://ec-jrc.github.io/lisflood-code/4_Static-Maps_soil-hydraulic-properties/.
2. Assigning 0 to Ksat1 value is not the correct solution. Table 6 https://hess.copernicus.org/articles/28/2991/2024/#section13&gid=1&pid=1 explains how to assign values when the fraction is zero. Specifically, please refer to the third column (No Data filling).

It is our pleasure to support your implementation of OS LISFLOOD,
we hope that our answers clarified your doubts,
kind regards,
Stefania

[Hrekhani97]

Dear @StefaniaGrimaldi
Thank you for the clarification provided earlier. I referred to the document and successfully created all the soil hydraulic parameters.

I am currently working on preparing the channel geometry layers and have encountered a few questions. In some areas, such as the downstream floodplain, the channel gradient value is 0, while in mountainous regions, the gradient is greater than 1. However, the website states that the channel gradient should fall within the [0,1] range. At the same time, the Channel Gradient or Slope map for the European domain (at 1 arc-minute horizontal resolution, shown on the left) has values ranging from 0 to 1.58109. This has left me a bit confused. Could you please clarify this discrepancy?

Thank you for your continuous support. I look forward to your guidance on this matter.

Kind regards

[StefaniaGrimaldi]

Dera @Hrekhani97,

congratulations on the successful preparation of the soil hydraulic parameters.

Could you please help me understand your recent enquiry about the channel geometry layers?

Your question concerns the channel gradient map explained in https://ec-jrc.github.io/lisflood-code/4_Static-Maps_channel-geometry/, correct?
I downloaded the 1 arcmin map changrad_European_01min.nc from https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/CEMS-EFAS/LISFLOOD_static_and_parameter_maps_for_EFAS/Catchments_morphology_and_river_network/ and I could not identify the pixels that you mentioned.
In your message, I read "shown on the left": could you please clarify where should I look?

Regarding the 0 values, please note that the LISFLOOD .xml settings file includes the key "ChanGradMin": https://github.com/ec-jrc/lisflood-code/blob/master/src/lisfloodSettings_reference.xml#L1442
You can see the usage of this latter key within the kinematic wave computations here: https://github.com/ec-jrc/lisflood-code/blob/master/src/lisflood/hydrological_modules/routing.py#L184

Kind regards,
Stefania