srlabUsask / crhmcode

GNU General Public License v3.0
8 stars 4 forks source link

glacier module debris-cover melt issue #432

Closed loganxingfang closed 5 months ago

loganxingfang commented 1 year ago

I noticed the lagT_used and lagSW_used are not simulated correctly in glacier module in CRHMcode version 1.3.3. I attached a prj file along with the obs file to demonstrate the issue.

The prj file is here PeytoDebrisCoverTest_glacier5_newObs.zip The obs file is here ERAI_teauQsiQlip_MonthlyAdjust2peytoBow_Jul2017_DEB_newDateFormat.zip

  1. The simulated lagT_used and lagSW_used are different between CRHM_GUI and Borland CRHM:

Screen shots for these simulated by CRHM_GUI: lagT_used_CRHM_GUI lagSW_used_CRHM_GUI

Screen shots for these simulated by Borland CRHM: lagT_used_BorlandCRHM lagSW_used_BorlandCRHM

  1. In addition, I output the lagT_delayed and lagSW_delayed from both CRHM_GUI and Borland CRHM, as they are directly assigned to lagT_used and lagSW_used. I only output the first layer for all 7 HRUs for demonstration, and there are 11 other layers that cannot be displayed in one GUI. Anyway, the simulated lagT_delayed and lagSW_delayed are not correct.

Screen shots for these simulated by CRHM_GUI: lagT_delayed_CRHM_GUI lagSW_delayed_CRHM_GUI

Screen shots for these simulated by Borland CRHM: lagT_delayed_BorlandCRHM lagSW_delayed_BorlandCRHM

  1. Next, the layer index are simulated as lagT and lagSW based on debris_h and use_debris parameters. I output lagT and lagSW from both CRHM_GUI and Borland CRHM. Both lagT and lagSW are correctly simulated in CRHM_GUI.

Screen shots for these simulated by CRHM_GUI: lagT_CRHM_GUI lagSW_CRHM_GUI

Screen shots for these simulated by Borland CRHM: lagT_BorlandCRHM lagSW_BorlandCRHM

  1. Furthermore, I noticed 6/1/2017 is shown as start date in the lower left TChart Panel of CRHM_GUI, while it should be just a few hours prior to 6/30/2017 as shown in Borland CRHM. I checked the CRHM_output from CRHM_GUI, and in that the first time step is correct: 6/30/2017 1:00:00 AM. I am not sure why the display is not correct in the TChart Panel, and I have screen shots below:

start date in lower left corner of TChart Panel of CRHM_GUI: start_date_CRHM_GUI

start date in lower left corner of TChart Panel of Borland CRHM: start_date_BorlandCRHM

jhs507 commented 5 months ago

Could you upload the outputs of these variables as produced by BorlandCRHM I don't have it on my machine. I think looking at the specific values will be helpful here.

loganxingfang commented 5 months ago

Could you upload the outputs of these variables as produced by BorlandCRHM I don't have it on my machine. I think looking at the specific values will be helpful here.

I attached the Borland CRHM output file here CRHM_output_1.txt

jhs507 commented 5 months ago

I am still investigating but in Borland CRHM LagT and LagSW are floating point variables and in crhmcode they are both integer variables.

This could be the root of the issue should they be floating point variables or integer variables.

jhs507 commented 5 months ago

Also in both crhmcode and borland the calculation calls a function called depth of snow. In both borland and crhmcode there is a local definition of depth of snow in the classglacier class but in crhmcode the depth of snow function being called is the one found in common instead.

In borland and crhmcode the depth of snow function in the glacier module is:

float Classglacier::DepthofSnow(float SWE){ // (mm)

// Tabler et al. (1990b) Calculates Snow Depth(mm) from SWE(mm)
  float rho;

  if (SWE > 1.0) {
    rho = 522.0 - 204700/SWE*(1.0 - exp(-SWE/673.0)); // converted from original cm to mm
    return SWE/rho*1000.0; // (mm)
  }
  else
    return 0.0;

} // DepthofSnow

In crhmcode the common depth of snow function is:

double Common::DepthofSnow(double SWE)
{
    /* 3/5/98
    Calculates Snow Depth(m) from SWE(mm) */
    double Snow_Depth;

    if (SWE > 2.05) {
        if (SWE <= 145.45) /* SWE 145.45 mm equivalent to 60 cm*/
            Snow_Depth = (SWE - 2.05) / 2.39;
        else
            Snow_Depth = (SWE + 128.06) / 4.5608;
    }
    else
        Snow_Depth = 0;

    return Snow_Depth / 100.0;
} /* DepthofSnow*/SS

The common DepthofSnow is being used in the glacier module in crhmcode.

loganxingfang commented 5 months ago

Also in both crhmcode and borland the calculation calls a function called depth of snow. In both borland and crhmcode there is a local definition of depth of snow in the classglacier class but in crhmcode the depth of snow function being called is the one found in common instead.

In borland and crhmcode the depth of snow function in the glacier module is:

float Classglacier::DepthofSnow(float SWE){ // (mm)

// Tabler et al. (1990b) Calculates Snow Depth(mm) from SWE(mm)
  float rho;

  if (SWE > 1.0) {
    rho = 522.0 - 204700/SWE*(1.0 - exp(-SWE/673.0)); // converted from original cm to mm
    return SWE/rho*1000.0; // (mm)
  }
  else
    return 0.0;

} // DepthofSnow

In crhmcode the common depth of snow function is:

double Common::DepthofSnow(double SWE)
{
  /* 3/5/98
  Calculates Snow Depth(m) from SWE(mm) */
  double Snow_Depth;

  if (SWE > 2.05) {
      if (SWE <= 145.45) /* SWE 145.45 mm equivalent to 60 cm*/
          Snow_Depth = (SWE - 2.05) / 2.39;
      else
          Snow_Depth = (SWE + 128.06) / 4.5608;
  }
  else
      Snow_Depth = 0;

  return Snow_Depth / 100.0;
} /* DepthofSnow*/SS

The common DepthofSnow is being used in the glacier module in crhmcode.

In Classglacier.cpp file, there is code calculating DepthofSnow from line 962 to 974: double Classglacier::DepthofSnow(double SWE) { // (mm)

// Tabler et al. (1990b) Calculates Snow Depth(mm) from SWE(mm) double rho;

if (SWE > 1.0) {
    rho = 522.0 - 204700 / SWE * (1.0 - exp(-SWE / 673.0)); // converted from original cm to mm
    return SWE / rho * 1000.0; // (mm)
}
else
    return 0.0;

} // DepthofSnow

The above is same as that in Borland CRHM and should be called directly by Classglacier in crhmcode; this can be done by changing the code on line 881 to read this subroutine instead of that from Common.

glacier_depth[hh] = (Classglacier::DepthofSnow(SWE[hh]) + firn_depth[hh]) / 1000.0 + ice[hh] / ice_dens[hh]; // (m) thickness

jhs507 commented 5 months ago

What about LagT and LagSW being integer values instead of floating point values? Is this an issue as well?

loganxingfang commented 5 months ago

What about LagT and LagSW being integer values instead of floating point values? Is this an issue as well?

Both lagT and lagSW should be declared as long, which is correctly done in crhmcode. Their values reflecting hours lagged behind inputs of T and SW as result of debris cover above glacier.

jhs507 commented 5 months ago

I have switched the glacier module to use the locally declared version of DepthofSnow but I still get the same values for the variables as before. And do not get the oscillating behavior in your borland CRHM screenshots.

loganxingfang commented 5 months ago

I checked the code again, and there is nothing really different in crhmcode compared to that in Borland CRHM.

The point of this debris cover glacier to lag T and SW. The steps are:

First using the two-dimensional variables to store value by moving data up, or reverse order; this should be achieved by the following code: from line 532 to 536 // 22 Apr 2022 modified: for(long ll = NUM_SIZE - 1; 0 < ll; --ll){ // move data up lagT_delayed_lay[ll][hh] = lagT_delayed_lay[ll-1][hh]; lagSW_delayed_lay[ll][hh] = lagSW_delayed_lay[ll-1][hh]; }

Then, get new value for these: from line 537 to 539 // enter new data lagT_delayed_lay[0][hh] = hru_t[hh]; lagSW_delayed_lay[0][hh] = Qnsn_Var[hh];

Next, based on value of lagT and lagSW as result of debris cover, assign values of lagT_delayed_lay and lagSW_delayed_lay and use them for calculating enery flux: from line 540 to 544 // lag T and SW long LT = lagT[hh]; long LSW = lagSW[hh]; lagT_used[hh] = lagT_delayed_lay[LT][hh]; lagSW_used[hh] = lagSW_delayed_lay[LSW][hh];

Based on the screen shots I included for the issue, crhmcode appears unable to generate correct values of lagT_used and lagSW_used. Maybe there is something preventing correct LT and LSW being used for lagT_delayed_lay[LT][hh] and lagSW_delayed_lay[LSW][hh].

jhs507 commented 5 months ago

I have taken a closer look and determined that something else is likely going on.

I don't know what the cause is but instead of stepping through hourly from 6/30/2017 to 7/08/2017 only two time steps are output and then the model stops running leading to the leading to outputs being entirely wrong.

6/1/2017 1:00:00 AM 7/4/2017 1:00:00 AM

Are the two time steps output. I don't know the cause of this behavior but it will certainly affect more than the variables we are investigating.

jhs507 commented 5 months ago

Found a possible fix and the results look promising.

2024-05-28 15_55_34-The Cold Regions Hydrological Model - C__Users_jhs507_repos_crhmcode_crhmcode_pr

jhs507 commented 5 months ago

Okay I think I have found the issue. If the start date of the simulation and the date of the first observation are not aligned a one line bug was causing the issues I discovered earlier with only two time steps to be output. Fixing this bug may or may not fix glaicer debris-cover melt issue as we have not actually been looking a proper data for it.

I am going open a new issue for this bug and fix it.

loganxingfang commented 5 months ago

Okay I think I have found the issue. If the start date of the simulation and the date of the first observation are not aligned a one line bug was causing the issues I discovered earlier with only two time steps to be output. Fixing this bug may or may not fix glaicer debris-cover melt issue as we have not actually been looking a proper data for it.

I am going open a new issue for this bug and fix it.

@jhs507 Yes, the start and end date for model simulation is defined in the prj file under Dates:

Dates:

2017 6 30 2017 7 8

These above are defined by user either entering it outside CRHM GUI using text edit or using select button from Start Date and End Date.

The start and end date are not always same and defined by the date length in observation data file.

jhs507 commented 5 months ago

Yeah, the bug is described here #448

jhs507 commented 5 months ago

After the DepthofSnow fix and fixing issue 448 These are the results I get.

TestOutput.obs.txt TestOutput.txt

jhs507 commented 5 months ago

These results seem to be comparable to me. I will compile a version to send to you for testing.

jhs507 commented 5 months ago

@loganxingfang You can use this version to test to see if it solves this issue.

Install_CRHM_GUI_issue_432.zip

loganxingfang commented 5 months ago

@loganxingfang You can use this version to test to see if it solves this issue.

Install_CRHM_GUI_issue_432.zip

@jhs507 The output file CRHM_output_1.zip from CRHM_GUI_issue_432 version is almost same as that CRHM_output_1.txt from Borland CRHM. There are some noticeable different values between them on or around July 1.

loganxingfang commented 5 months ago

@loganxingfang You can use this version to test to see if it solves this issue. Install_CRHM_GUI_issue_432.zip

@jhs507 The output file CRHM_output_1.zip from CRHM_GUI_issue_432 version is almost same as that CRHM_output_1.txt from Borland CRHM. There are some noticeable different values between them on or around July 1.

@jhs507 I found a bug in glacier module debris-cover melt, this is likely the cause for the difference. See below:

lagSW_delayed_lay[0][hh] = Qnsn_Var[hh];

The above should be:

lagSW_delayed_lay[0][hh] = Qsisn_Var[hh];

It is the incoming shortwave radiation variable Qsisn_Var that is being delayed. Not the net radiation variable Qnsn_Var. I will make correction in Borland CRHM and also create a branch for srlab crhmcode to be pulled.

jhs507 commented 5 months ago

I already have a branch for this fix you can apply the borland changes to this branch and we can merge it into master as one pull request. https://github.com/srlabUsask/crhmcode/tree/issue_432_glacier_debris_cover_melt

loganxingfang commented 5 months ago

I already have a branch for this fix you can apply the borland changes to this branch and we can merge it into master as one pull request. https://github.com/srlabUsask/crhmcode/tree/issue_432_glacier_debris_cover_melt

Will do.

jhs507 commented 5 months ago

I have made the change in CRHMcode dfdda8741be03be3374edf6e6e777e150e202a75 with this commit

loganxingfang commented 5 months ago

I have made the change in CRHMcode dfdda87 with this commit

Great. I uploaded updated Borland source code files to borland_source folder on the main branch.

jhs507 commented 5 months ago

With these fixes is this issue fixed?

loganxingfang commented 5 months ago

With these fixes is this issue fixed?

In Borland CRHM, yes. I will take a look at you CRHM_GUI if you send me an installation file.

loganxingfang commented 5 months ago

@loganxingfang You can use this version to test to see if it solves this issue.

Install_CRHM_GUI_issue_432.zip

@jhs507 I also used this version for running a researcher's prj files. I have some issues. Should I open an issue for that?

jhs507 commented 5 months ago

If the issues are not related to this issue then opening a new issue would be good.

jhs507 commented 5 months ago

Install_CRHM_GUI_issue_432_v2.zip

Try this to see if it solves the debris-cover melt issues. If you have discovered other issues it would be helpful to open a new issue unless the issues are specifically caused by the fix above and don't occur in the master branch.

loganxingfang commented 5 months ago

Install_CRHM_GUI_issue_432_v2.zip

Try this to see if it solves the debris-cover melt issues. If you have discovered other issues it would be helpful to open a new issue unless the issues are specifically caused by the fix above and don't occur in the master branch.

Yes, this is fixed. CRHM_GUI_issue_432_v2 resolved the issue.

jhs507 commented 5 months ago

Alright I will merge this in and then be releasing a new version with the latest fixes.