0

possibly bug in viscosity calculation in TOUGH3_EWASG

Dear developers, 

 

First of all, I would like to say thank you for developing this excellent codes.

Here I am reporting a tiny bug I found in viscosity calculation in TOUGH3 EWASG, as I could not find any relevant report in other posts. I'd be sorry if here is not right place to report a bug, or the issue is already solved.

 

When I set IE(3)=0 and IE(4)=6 (Driesner 2007), TOUGH3 calculated viscosity values are very small (order of 1E-5) while fluid is in liquid phase and expected value is order of 1E-4. 

After quickly checking the code, I found that fluid density (DW0) is probably not correctly updated before calculating water viscosity if IE(4)=6. I guess the density variable remains its initial value of zero and this causes the low viscosity.

 

Kind regards,

Norihiro

6 replies

null
    • kenny
    • 2 wk ago
    • Reported - view

    Hi Norihiro,

    Thanks for reporting the potential bug. I am very interested in investigating this problem. If you can provide your input files, it will be a great help. 

    Kenny

      • Norihiro_Watanabe
      • 2 wk ago
      • Reported - view

       

      Thank you for the reply. I'm very glad that you are interested in this potential bug.  

      I can provide you two input files (attached). One is with default case (IE(4)=0) and the other one is with the Driesner model (IE(4)=6). Problem settings in the files are based on the 1D salt water multiphase simulation in the following paper (I was trying to reproduce their Fig 6 although it was not successful with configuration written in the paper). 

      Weis, P., Driesner, T., Coumou, D., & Geiger, S. (2014). Hydrothermal, multiphase convection of H2O-NaCl fluids from ambient to magmatic temperatures: a new numerical scheme and benchmarks for code comparison. Geofluids, 14(3), 347–371.

      Running TOUGH3 with the two input files, you can clearly see difference in calculated liquid viscosity (VIS_L) in CSV output.

       

      Best regards,

      Norihiro 

      INFILE_defaultINFILE_driesner

    • kenny
    • 2 wk ago
    • Reported - view

    Hi Norihiro,

    You are absolutely correct. The problem is caused by not updating the DW0 for the case IE(4)=6.

    This problem can be fixed by insert following two lines into (bold lines) the subroutine  COBRI in file thrmb.f90 (at around line 290):

          IF(IE(4).EQ.6) THEN
            CALL DRIESNER(TX,PX,XS,DB,HB)
            PL=MAX(PX,PSATW)
            CALL COWAT(TX,PL,DW0,UB)         
            RETURN
          ENDIF

    Thanks for your contributions making TOUGH more robust.

     

    Kenny

      • Norihiro_Watanabe
      • 2 wk ago
      • Reported - view

       Thank for investigating the issue and posting the fix. 

      When I tried to fix the code following your suggestion, I got one question. In COBRI() in thrmb.f90, it looks that, if IE(4) is not 6, DW0 is set as a liquid density at water saturation pressure. Below is example code in the subroutine.

      CALL COWAT(TX,PSATW+DP,DW,UW)

      DW0=DW

      On the other hand, in your above post, DW0 is calculated from original pressure, not saturation pressure. I wonder which one is correct, though it might depend on a viscosity function to be used. Or is it assumed that the pressure difference is neglegible for liquid density?

      • kenny
      • 2 wk ago
      • Reported - view

       My understanding is that the pressure used for the liquid density calculation is liquid water pressure. COBRI() was designed for both single aqueous phase or two-phase flow. In single phase condition, PX must be used, in two-phase condition the saturation pressure must be used. 

      • Norihiro_Watanabe
      • 2 wk ago
      • Reported - view

       Thank you for the clear answer. I agree with you  

Login to reply

Content aside

  • 2 wk agoLast active
  • 6Replies
  • 18Views
  • 2 Following