New Green Ampt Method

[FLO-2DKaren]

The DTHETA calc on the new Green Ampt Method is wrong. We need to adjust the landuse soil intersection so that...

Intersect SOIL to LandUse and calculate get dTHETA parts.

if ISAT = wet, DTHETA = 0, if ISAT = normal DTHETA part is Normal if ISAT = dry DTHETA part is Dry.

Here is a sample of the data. It goes along with the Land use from Lesson 1 Hydrology. AZ Green Ampt.zip

[FLO-2DKaren]

@rpachaly Hi Robson, I updated this job. I just realized we made a mistake on the current infiltration editor. Let's have a meeting when you can and I'll show you the issue.

[FLO-2DKaren]

The select fields form will need DTHETAdry and DHTEDAnormal fields selections.

[rpachaly]

Hi Karen,

Anytime after 1230 your time works for me.

[FLO-2DKaren]

The code is in infiltration_tools.py near line 164 You can see that the land_soil_geom is the location where the new field will be added.

for land_feat, engine in land_features.values():
                    land_rtimp = land_feat[self.rtimpl_fld]
                    land_saturation = land_feat[self.saturation_fld]
                    land_geom = land_feat.geometry()

                    soil_fids = soil_index.intersects(land_geom.boundingBox())
                    for soil_fid in soil_fids:
                        soil_feat, soil_engine = soil_features[soil_fid]
                        soil_rtimp = soil_feat[self.rtimps_fld]
                        soil_dtheta = soil_feat[self.soil_dtheta_fld]
                        soil_psif = soil_feat[self.soil_psif_fld]

                        land_soil_geom = land_geom.intersection(soil_feat.geometry())

                        if land_soil_geom.isEmpty():
                            continue

                        land_soil_feat = QgsFeature(land_soil_fields, land_soil_fid)

                        land_soil_feat.setGeometry(land_soil_geom)
                        land_soil_feat["rtimp"] = max(land_rtimp, soil_rtimp)
                        land_soil_feat["dtheta"] = soil_dtheta
                        land_soil_feat["psif"] = soil_psif
                        land_soil_feat["saturation"] = land_saturation
                        land_soil_features[land_soil_fid] = (
                            land_soil_feat,
                            QgsGeometry.createGeometryEngine(land_soil_geom.constGet()),
                        )
                        land_soil_index.addFeature(land_soil_feat)
                        land_soil_fid = land_soil_fid + 1

[FLO-2DKaren]

The dialog box code is in infil_green_ampt.ui. The dialog connection code is in infil_editor_widget line 789

uiDialog_green, qtBaseClass_green = load_ui("infil_green_ampt")

The code for the processing is in infiltration_tools.py

[FLO-2DKaren]

AZ Green Ampt 2223.zip

For Testing... Use Lesson 1 Part 1 to get the land use data. Set several of the polygons to "wet" Use the shapefile attached to this comment for the soil data.

The test methods video is here: https://transcripts.gotomeeting.com/#/s/5e14b8a7e6aeb513119d962a44e0910471ef6ea9252624a458c507172562997b

[rpachaly]

@FLO-2DKaren Hi Karen,

I'm testing my modifications on the code here and one question arose. If I don't check the the Log Average Area, does it use the dtheta (normal or dry) in any calculations? As far as I understood the code, it uses only the avg_xksat and saturation to determine the dtheta using the calculate_dtheta method.

[FLO-2DKaren]

@rpachaly Yes that is correct. That is the 2018 method. That method is OK.

We need to modify the 2023 method only. The other method is OK.

[rpachaly]

Ok, good!

[rpachaly]

@FLO-2DKaren

I finished the New Green-Ampt method as requested. Everything looks alright. My testing procedure was the following. I used the lesson 1 project with the land use and soil data provided. I changed two land use areas to “wet” and “saturated”. Then, I joined the “Grid” shapefile (under the Schematic Layers) with “Cells Green-Ampt” to see if each cell was returning the expected value. I calculated the Green-Ampt three times:

• Setting dtheta normal and dtheta dry on the Calculate Green-Ampt to the field DTHETANrm. In this case, it should return only the dtheta normal value. Which returned the expected value. • Setting dtheta normal and dtheta dry on the Calculate Green-Ampt to the field DTHETADry. In this case, it should return only the dtheta dry value. Which returned the expected value. • Using the correct fields for dtheta normal and dtheta dry on the Calculate Green-Ampt. In this case, it should return the log area average. Which returned the expected value.

Finally, I checked if all “wet” and “saturated” areas were returning 0 as requested, which was also correct.

However, I have two questions here:

1. Since we’re using a Log Area Average, when we have a cell with very high wet/saturated area, the dtheta will be close to 1. Then, when the cell is entirely wet/saturated, the dtheta is 0. Is that correct? Here is an example: Let’s suppose that 99% of a cell is wet (dtheta = 0) and 1% is normal (dtheta = 0.25). Let’s assume an area of 1 m2.

   Avg_dtheta = 0.99*0 + 0.01*log(0.25)
   Avg_dtheta = -0.006
   Avg_dtheta = 10^(-0.006)
   Avg_dtheta = 0.98

2. If the previous question is true, in some cells, the dtheta will be 0.999 and it could be automatically rounded to 1. Is that an issue?

If these questions are not relevant, I can push the branch.

In the meanwhile, I’ve receiving some DeprecationWarnings here. I’ll solve them.

[FLO-2DKaren]

@rpachaly That is certainly a relevant issue. It means we are applying the dtheta = zero incorrectly for any cell that is partly wet and partly dry or normal.

What are your thoughts on this issue. I know there's a simple solution but I can't remember how we did it in the past.

[rpachaly]

Well, I didn't quite understand the reason for using a log average area in this case. In my opinion, I would just use a weighted average for computing the average dtheta in a cell. In the code, we have the area and the dtheta inside the for loop. It is really simple to implement the weighted average.

[rpachaly]

Let me know if I can proceed with this approach 🙂

[FLO-2DKaren]

@rpachaly @vaughanster

Haha Yes this is actually a very common question. So Nate told me that the xksat of the soil layer statistically is a log distribution. Test results and analyses show that it is the more accurate method. They have a writeup in their hydrology documentation that shows how this works.

I'm more a fan of not averaging it at all and letting the grid size be the control for the sample. FCDMC won't have it!!! So here we are.

[FLO-2DKaren]

Nate is going to show us how he does this using ArcGIS.

[FLO-2DKaren]

@rpachaly Hi Robson,

I found this in a email from Nate. It looks like we simply leave the 0 part out of the equation for any of the psif, xksat, or dtheta. I don't think I have ever seen those values except for dtheta. I can look at my mapping data and see if it exists.

if DTHETA = 0 (saturated soil), the infiltration rate becomes mathematically equal to XKSAT. While PSIF isn't 0 in that case, it might as well be. Calculating the infiltration depth become impossible in this case as well (which Jim and Noemi are well aware).

What you describe between the soil and land-use layers is the intended use - the potential DTHETA values are defined by the soil, but the operational antecedent moisture condition is defined by the land-use.

In general, the XKSAT log-averaging method should be okay for PSIF and DTHETA as well - handling of 0 values, which you've already touched on, is the key difference. The easy way to handle the 0 values for DTHETA is to exclude them from the log-average numerator, but not the area from the denominator

That's really what the "if xksat > 0" part at the end of this line is doing: xksat_gen = [area * log10(xksat) for xksat, area in parts if xksat > 0]

Because sub-part area is calculated in %, the reduction is handled automatically by omitting any sub-parts with an xksat (or DTHETA, PSIF, etc.) with a value of 0. It's pretty elegant, really. Care should be taken to handle the case where xksat_gen has no non-zero sub-parts (e.g. the whole thing has a 0 value) because the statement above isn't evaluatable in that case - I think xksat_gen would resolve to None (null) in that case and would need to be set to 0.

[vaughanster]

Karen,

I don't know that this holds up mathematically for DTHETA, which can reasonable be 0. If XKSAT goes to 0, there is no infiltration and this should be modeled as RTIMP - we simply don't calculate XKSAT or PSIF == 0 in a normal condition. DTHETA, by comparison, is reasonably 0 - my thinking currently is that DTHETA == 0 components should be excluded from the logarithmic areal weighting and calculated separately using simple areal weighting with the non-zero components calculated with logarithmic weighting as directed in the FCDMC manual.

Nate Vaughan, PE JE Fuller Hydrology & Geomorphology, Inc. 143 N McCormick Street, Suite 203 Prescott, AZ 86301

@.*** 928-640-0778

Engineering data, calculations, or opinion transmitted via e-mail are Preliminary and Not For Construction.

---

From: Karen @.> Sent: Saturday, June 24, 2023 6:29 AM To: FLO-2DSoftware/qgis-flo-2d-plugin @.> Cc: Nate Vaughan @.>; Mention @.> Subject: Re: [FLO-2DSoftware/qgis-flo-2d-plugin] New Green Ampt Method (Issue #959)

CAUTION: [EXTERNAL] : this e-mail originated from outside your organization. Exercise caution when opening attachments or clicking links, especially from unknown senders.

@rpachalyhttps://github.com/rpachaly Hi Robson,

I found this in a email from Nate. It looks like we simply leave the 0 part out of the equation for any of the psif, xksat, or dtheta. I don't think I have ever seen those values except for dtheta. I can look at my mapping data and see if it exists.

if DTHETA = 0 (saturated soil), the infiltration rate becomes mathematically equal to XKSAT. While PSIF isn't 0 in that case, it might as well be. Calculating the infiltration depth become impossible in this case as well (which Jim and Noemi are well aware).

What you describe between the soil and land-use layers is the intended use - the potential DTHETA values are defined by the soil, but the operational antecedent moisture condition is defined by the land-use.

In general, the XKSAT log-averaging method should be okay for PSIF and DTHETA as well - handling of 0 values, which you've already touched on, is the key difference. The easy way to handle the 0 values for DTHETA is to exclude them from the log-average numerator, but not the area from the denominator

That's really what the "if xksat > 0" part at the end of this line is doing: xksat_gen = [area * log10(xksat) for xksat, area in parts if xksat > 0]

Because sub-part area is calculated in %, the reduction is handled automatically by omitting any sub-parts with an xksat (or DTHETA, PSIF, etc.) with a value of 0. It's pretty elegant, really. Care should be taken to handle the case where xksat_gen has no non-zero sub-parts (e.g. the whole thing has a 0 value) because the statement above isn't evaluatable in that case - I think xksat_gen would resolve to None (null) in that case and would need to be set to 0.

— Reply to this email directly, view it on GitHubhttps://github.com/FLO-2DSoftware/qgis-flo-2d-plugin/issues/959#issuecomment-1605499857, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEHRXKEHF3K4ZS2MJO6SKJTXM3TUFANCNFSM6AAAAAAZNI7G5A. You are receiving this because you were mentioned.Message ID: @.***>

---

From: Karen @.> Sent: Saturday, June 24, 2023 6:29 AM To: FLO-2DSoftware/qgis-flo-2d-plugin @.> Cc: Nate Vaughan @.>; Mention @.> Subject: Re: [FLO-2DSoftware/qgis-flo-2d-plugin] New Green Ampt Method (Issue #959)

CAUTION: [EXTERNAL] : this e-mail originated from outside your organization. Exercise caution when opening attachments or clicking links, especially from unknown senders.

@rpachalyhttps://github.com/rpachaly Hi Robson,

I found this in a email from Nate. It looks like we simply leave the 0 part out of the equation for any of the psif, xksat, or dtheta. I don't think I have ever seen those values except for dtheta. I can look at my mapping data and see if it exists.

if DTHETA = 0 (saturated soil), the infiltration rate becomes mathematically equal to XKSAT. While PSIF isn't 0 in that case, it might as well be. Calculating the infiltration depth become impossible in this case as well (which Jim and Noemi are well aware).

What you describe between the soil and land-use layers is the intended use - the potential DTHETA values are defined by the soil, but the operational antecedent moisture condition is defined by the land-use.

In general, the XKSAT log-averaging method should be okay for PSIF and DTHETA as well - handling of 0 values, which you've already touched on, is the key difference. The easy way to handle the 0 values for DTHETA is to exclude them from the log-average numerator, but not the area from the denominator

That's really what the "if xksat > 0" part at the end of this line is doing: xksat_gen = [area * log10(xksat) for xksat, area in parts if xksat > 0]

Because sub-part area is calculated in %, the reduction is handled automatically by omitting any sub-parts with an xksat (or DTHETA, PSIF, etc.) with a value of 0. It's pretty elegant, really. Care should be taken to handle the case where xksat_gen has no non-zero sub-parts (e.g. the whole thing has a 0 value) because the statement above isn't evaluatable in that case - I think xksat_gen would resolve to None (null) in that case and would need to be set to 0.

— Reply to this email directly, view it on GitHubhttps://github.com/FLO-2DSoftware/qgis-flo-2d-plugin/issues/959#issuecomment-1605499857, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEHRXKEHF3K4ZS2MJO6SKJTXM3TUFANCNFSM6AAAAAAZNI7G5A. You are receiving this because you were mentioned.Message ID: @.***>

[rpachaly]

I don't know that this holds up mathematically for DTHETA, which can reasonable be 0. If XKSAT goes to 0, there is no infiltration and this should be modeled as RTIMP - we simply don't calculate XKSAT or PSIF == 0 in a normal condition. DTHETA, by comparison, is reasonably 0 - my thinking currently is that DTHETA == 0 components should be excluded from the logarithmic areal weighting and calculated separately using simple areal weighting with the non-zero components calculated with logarithmic weighting as directed in the FCDMC manual.

How can I have access to the FCDMC manual? I want to take a look to fully understand what is being discussed here.

[FLO-2DKaren]

Hi Robson,

I'm OK with your method for DTHETA = 0 arithmetic average. The method should be described in chapter 4 but it might be for HEC-1 not for FLO-2D.

Chapter 4 - Rainfall Losses.pdf

[FLO-2DKaren]

@rpachaly Robson,

I think they are waiting for us to figure this out so they can adapt their documentation to our method. So we'll need approval for anything we do but we should have our own ability to make the decision.

[rpachaly]

Okay, let me read the key points of this chapter and think a little bit on how to proceed.

[FLO-2DKaren]

I also think Nate's last email is relevant too.

[rpachaly]

My opinion: I agree with Nate. We need to use logarithmic areal weighting for dry & normal areas in a cell. When we have wet, dry and normal in a cell, we need to calculate separately the wet portion using simple areal weighting. In this case, we will have two steps for calculating the dtheta.

1. Calculate the equation 4.4 on chapter 4 considering ONLY the area not saturated (dry & normal)

3. If no wet areas are present, dtheta_dn is dtheta_avg. If there are wet areas, use the following equation:

I ran some scenarios on excel to see if everything holds up and I believe that this approach works for dtheta. Let me know if you agree with this approach.

[FLO-2DKaren]

@vaughanster Nate, What do you think of Robson's approach?
I think I sent my email to Elizabeth too soon.

[FLO-2DKaren]

@rpachaly Hi Robson, I have some feedback from the Flood Control District regarding DTHETA. They are going to reassess their method and get back to us on that.

We don't have to wait for them because we can get our product out the door now and redo the code when they make a final decision.

For now, please change DTHETA method to an arithmetic mean so that the 0 parts don't mess up the equation.

[rpachaly]

Hi Karen,

Do you want me to get back on this issue or keep moving forward with the Infiltration Editor? I remember that you said something about a deadline on this Friday.

[FLO-2DKaren]

@rpachaly Hi Robson, Do you get emails when I don't "mention" you in the comment?

Let's finish the calculator first so we can publish something. Then we'll do the other stuff.

[rpachaly]

Hi Robson, Do you get emails when I don't "mention" you in the comment?

I get all emails, but I think only when I had been mentioned on the issue.

The calculator you mean this New Green Ampt Method?

[FLO-2DKaren]

@rpachaly Yes I meant to work on the DTHETA Calculator. Let's apply the arithmetic weighted avg first since it is so close to being ready. We can come back and modify it if the District wants something different.

[rpachaly]

I finished up the the DTHETA calculator. I've the two methods here ready to go, both tested and fully functional. One using only arithmetic weighted avg for Dry, Normal, and Wet saturations and the other using log area average for Dry and Normal and weighted average for Wet.

Which one do you want me to push the branch? I can also comment out one method and you select the one that you prefer.

[rpachaly]

I'll push the branch with the arithmetic weighted area for all saturations and the other method commented out in case we need to use something similar in the future.

[FLO-2DKaren]

@rpachaly Leave Both methods in the same code. Comment out the Log Weighted Average. We may adjust that in the future.

[FLO-2DKaren]

@rpachaly I'll run this test on UEFCC which is a huge model I have with Land Use data. I'll screenshot the test.

[rpachaly]

Okay, let me know if everything worked as expected. My tests were performed on the Lesson 1 project.

[FLO-2DKaren]

@rpachaly My calculator test is complete. It took 1 hour to process over 4 million cells. I think that is good.

Can you make a document that shows the validates the calculators on the small test case. We just need to validate the calculator for a variety of cells. change a few "wet" initial saturation conditions to the land use fields and then spot check the calculators. I think you can pick cells that have obvious blocks and approximate the areas. So long as we are within 0.01 on the calculator, we should be OK.

[rpachaly]

Ok! I'll do that.

[FLO-2DKaren]

@rpachaly Thanks Robson. That will help me a lot. I'm trying to wrap up another project right now.

[rpachaly]

G&A calculator validation.docx

I prepared this document. Let me know if this is what you were expecting.

I have to leave right now, but I'll be back in 2 hours. If you have any comments or suggestions, I can adjust the document then.

[FLO-2DKaren]

@rpachaly Yes that is perfect. Can you add the other calculators too?

xksat and psif are the same calculator so just test one of them.

rtimp is complex because it also does an intersection between the landuse and the soil before it intersects to the grid.

[rpachaly]

Hi Karen,

I didn't change the other calculators. PSIF uses the log area average. The RTIMP also uses the log average but it has some conditions for smaller areas.

Do you want me to change these calculators to weighted area average, or just add the other calculators as they are to the document?

[FLO-2DKaren]

@rpachaly I need the proof that we tested them for the District. Its one of our tasks.

[FLO-2DKaren]

@rpachaly I didn't answer your question. The log average calculators stay the same. Just need to test that they are functioning correctly.

[rpachaly]

Cool! I'll have that early this afternoon.

[rpachaly]

Hi Karen,

Let me know if I need to adjust, add, or remove anything.

G&A calculator validation.docx

[FLO-2DKaren]

@rpachaly

Looking good. I added a bit of detail to your document with respect to rtimp. It also uses a landsoil intersection for rtimpland and soilrockout. It takes the max of those two before it intersects to the grid for the weighted avg. Can you add that bit to your testing? Look for a rockout that is greater than rtimp land. You might have to modify the soil table to get one.

G.A.calculator.validation K.docx

Would you please review my Technical Reference Manual on the Green Ampt section. It needs an image with the calculator and look over my writeup and see if anything is wrong or could be written better.

FLO-2D Plugin Technical Reference Manual 2022.docx

I'm working on the QGIS User Manual. I'll ask you to look that over too when I am finished.

[rpachaly]

Here is the G&A calculator validation updated.

G.A.calculator.validation.Rob.docx

[rpachaly]

Here is the Technical Reference Manual. I added the two images for the G&A calculator and small changes on the text. I used the Track Changes.

FLO-2D.Plugin.Technical.Reference.Manual.2022.Rob.docx

Let me know if I can help in anything else. I'll be waiting for the QGIS User Manual.

[FLO-2DKaren]

OK the plugin portion of this is complete and tested. I'm going to close this issue. Great work.