Output values of MBio and MDispersed (in .sro file) are a function of number of particles (NBR_PARTIC)

[RachaelDMueller]

The notebook showing this problem is here

Created a new branch to work on

git checkout -b nbr_partic

NBR_PARTIC becomes NewOrigin%NbrParticlesIteration

My understanding from conversations with Shihan was that I ought to ignore ModuleLagrangian.F90 and that only ModuleLagrangianGlobal.F90 is used; but I'm not seeing that as true in the code. First, both are compiled. Second, if I search NBR_PARTIC in the code, I only see it show up in ModuleLagrangianGlobal.F90 within a call for NewOrigin%State%Deposition and under a call for !Min Sedimentation velocity. We aren't including sedimentation, so this call doesn't seem to apply to our case. In contrast, within ModuleLagrangian.F90, I see:

      !Reads parameter specific to cada Spatial emission type
PA:   if (NewOrigin%EmissionSpatial == Point_ .or.                             &
          NewOrigin%EmissionSpatial == Accident_ ) then

          !Gets the number of particles to emit
          call GetData(NewOrigin%NbrParticlesIteration,                        &
                       Me%ObjEnterData,                             &
                       flag,                                                   &
                       SearchType   = FromBlock,                               &
                       keyword      ='NBR_PARTIC',                             &
                       ClientModule ='ModuleLagrangian',                       &
                       Default      = 1,                                       &
                       STAT         = STAT_CALL)
          if (STAT_CALL /= SUCCESS_) stop 'ConstructOrigins - ModuleLagrangian - ERR990'

Unfortunately, when I look up NewOrigin%NbrParticlesIteration, I run into a dead end as it appears to only affect the volume in the case of a box spill.

   BX:     if (NewOrigin%EmissionSpatial == Box_) then
   ...
   if (flag==0) NewOrigin%ParticleBoxVolume = NewOrigin%PointVolume / NewOrigin%NbrParticlesIteration

We are using accident. This is our setup:

NBR_PARTIC                : 2000
EMISSION_SPATIAL          : Accident
EMISSION_TEMPORAL         : Instantaneous
POINT_VOLUME              : 1000

Looking up EmmissionAccident, I see:

!Calcula a area ocupada por cada tracador
        AreaParticle = AreaTotal / float(CurrentOrigin%NbrParticlesIteration)
...
   NewParticle%Geometry%Volume         = CurrentOrigin%PointVolume /                &
                                              float(CurrentOrigin%NbrParticlesIteration)
   NewParticle%Geometry%VolVar         = 0.0
   !Stores initial Volume
   NewParticle%Geometry%InitialVolume = NewParticle%Geometry%Volume

NewOrigin and CurrentOrigin are just different pointer names passed into different subroutines. I’m not yet sure what the design motivation is in having different names for these pointers.

From this portion of the code, I see that the New Particle volume is calculated by dividing the total volume by the number of particles. Not seeing the problem here. Taking a different approach.

Evaluating evaporation parameterization
The mass in the Evaporation is consistent among the different NBR_PARTIC tests.
Looking into the difference between the Evaporation parameterization and Biodegredation.
The Evaporation method we used for these runs is Fingas. For our case, with Fingas evaporation using
coefficients, the equation is:

InternalConstant         = (Me%Var%MassINI/100.0)                           &
                                               * (Me%Var%Fingas_Evap_Const1                     &
                                               + Me%Var%Fingas_Evap_Const2                      &
                                               * Me%ExternalVar%WaterTemperature)

for logarithmic

Me%Var%MEvaporatedDT =  InternalConstant                                    &
                                                * exp(-Me%Var%MEvaporated/InternalConstant)/60.0

Same concept for square root. Note: use of Me%Var%MassINI.

[EDIT: This makes sense. MassINI is initial mass. MassOil is mass of oil left over after weathering. See below]

Evaluating dispersion parameterization
DISPERSIONMETHOD = NewDispersion

            Me%Var%MDispersedDT = (Me%Var%MassOil) / (1.0-Me%Var%MWaterContent)  &
                              *(1.0-Me%Var%VWatercontent) * P_Star *WCC           &
                               /(0.812*3.14*2.0*Me%ExternalVar%Wind/9.8)

In this case, Me%Var%MassOil is used rather than Me%Var%MassINI

These are set in OilPropIni in ModuleOil_0D.F90

 Me%Var%MassINI            = Me%Var%Density * Me%Var%VolInic
            Me%Var%MassOil            = Me%Var%MassINI
            Me%Var%VolumeOil          = Me%Var%VolInic
            do n=1,5
                Me%Var%SurfaceAnalyteMass(n)=Me%Var%MassOil*Me%Var%AnalytePercentage(n)*0.01
            end do

MassOil is set to MassINI

The difference is still noteworthy (and worth investigating the timing/call of OilPropIni) but I’m now wondering if this could be an output error. The two variables that are output are MBio and MEvaporated.

[RachaelDMueller]

Evaluate differences between MDispersed, Bio, and MEvaporated

Bio and MDispersed change with NBR_PARTIC; MEvaporated does not.

Values are written to .sro via the TimeSerieOutput subroutine in ModuleOil_0D.F90, lines 5311 - 5401. These are lines for MDispersed

Me%TimeSerie%DataLine(ColMEvaporated   ) = Me%Var%MEvaporated
Me%TimeSerie%DataLine(ColMDispersed    ) = Me%Var%MDispersed
Me%TimeSerie%DataLine(ColMBio          ) = Me%Var%MBio
...
call WriteTimeSerieLine(Me%ObjTimeSerie,                                        &
                                DataLine            = Me%TimeSerie%DataLine,            &
                                ExternalCurrentTime = ExternalCurrentTime,              &
                                STAT                = STAT_CALL)
        if (STAT_CALL /= SUCCESS_)                                                      &
                call SetError(FATAL_, INTERNAL_,                                        &
                             "Subroutine TimeSerieOutput - ModuleOil_0D. ERR01") 

Nothing fishy here. Go fish!

MDispersed

ModuleLagrangian.F90

if (Me%Var%OilBiodegradation)   then
        Me%Var%MDispersed                 = Me%Var%MDispersed - Me%Var%MBioDT+ Me%Var%MDispersedDT     &
                                                * Me%Var%DTOilInternalProcesses
 else
     Me%Var%MDispersed                 = Me%Var%MDispersed + Me%Var%MDispersedDT     &
                                                * Me%Var%DTOilInternalProcesses
 end if

MEvaporated

ModuleOli_0D.F90: Subroutine Evaporated

              if (Me%State%FirstStepIP)   then
                    
                    Me%Var%MEvaporatedDT = InternalConstant                                     &
                                               * sqrt(Me%Var%DTOilInternalProcesses/60.)        &
                                               / Me%Var%DTOilInternalProcesses
                
                else

                    Me%Var%MEvaporatedDT = InternalConstant * InternalConstant                  &
                                               / (2*Me%Var%MEvaporated*60.)

                end if   
                
                Me%Var%VEvaporatedDT     = Me%Var%MEvaporatedDT / Me%Var%Density*   &
                                         ((1.0-Me%Var%VWatercontent)/(1.0-Me%Var%MWatercontent))

...
if  (Me%Var%MassOil - (Me%Var%MEvaporatedDT) * Me%Var%DTOilInternalProcesses      &
           .GT. 0.0)  then  
            
            Me%Var%MassOil            = Me%Var%MassOil - (Me%Var%MEvaporatedDT)          &
                                            * Me%Var%DTOilInternalProcesses
            Me%Var%VolumeOil           = ((Me%Var%MassOil/(1.0-Me%Var%MWaterContent) )/   &
                                         Me%Var%Density)*(1.0-Me%Var%VWaterContent)
            !Me%Var%VolumeOil          = max(0.0,Me%Var%VolumeOil - (Me%Var%VEvaporatedDT)        &
             !                               * Me%Var%DTOilInternalProcesses)
        
        else    cd4
            Me%Var%MEvaporatedDT      = Me%Var%MassOil / Me%Var%DTOilInternalProcesses
            Me%Var%VEvaporatedDT      = Me%Var%MEvaporatedDT / Me%Var%Density*    &
                                     ((1.0-Me%Var%VWatercontent)/(1.0-Me%Var%MWatercontent))
            Me%Var%MassOil            = 0.0
            Me%Var%VolumeOil          = 0.0
        end if  cd4
...
Me%Var%MEvaporated    = Me%Var%MEvaporated + Me%Var%MEvaporatedDT                &
                                    * Me%Var%DTOilInternalProcesses

Aside

line 11187 shows that OilMass is set on the first time step to:

CurrentPartic%OilMass          = Me%ExternalVar%MassINI / CurrentOrigin%NbrParticlesIteration

[RachaelDMueller]

MassOil vs. MassINI

Exploring the difference between these two as I think it’s the cause of the NBR_PARTIC error lays in the difference between these two.

ModeulOil_0D.F90
While not a big bug, and not one we would likely notice, MassOil is currently adjusted for all weathering processes except biodegredation. I am not seeing an adjustment for MBio, as i am for all other weathering processes. e.g.:

Me%Var%MassOil            = Me%Var%MassOil - (Me%Var%MEvaporatedDT)          &
                                            * Me%Var%DTOilInternalProcesses

Me%Var%MassOil                = Me%Var%MassOil - (Me%Var%MdispersedDT)      &
                                                * Me%Var%DTOilInternalProcesses

Me%Var%MassOil            = Me%Var%MassOil - (Me%Var%MSedimentedDT)         &
                                            * Me%Var%DTOilInternalProcesses

Me%Var%MassOil          = Me%Var%MassOil - (Me%Var%MdissolvedDT)            &
                                              * Me%Var%DTOilInternalProcesses

Me%Var%MassOil          = Me%Var%MassOil - (Me%Var%MChemDispersedDT)    &
                                              * Me%Var%DTOilInternalProcesses

Me%Var%MassOil         = Me%Var%MassOil - (Me%Var%MOilRecoveredDT)      &
                                             * Me%Var%DTOilInternalProcesses

Although it IS accounted for in removing from dispersed mass:

       Me%Var%MDispersed                 = Me%Var%MDispersed - Me%Var%MBioDT+ Me%Var%MDispersedDT     &
                                                * Me%Var%DTOilInternalProcesses

subroutine OilInternalProcesses

MassINI is passed in as a variable to subroutine OilInternalProcesses, MassOil is not.

MassINI is set to Me%Var%MassINI
(where is Me%Var% specified?)

[RachaelDMueller]

Continuing on yesterday’s exploration of the usage of MassOil vs MassINI in the dispersion calculation. Yesterday, I experimented with using MassINI in place of MassOil in calculating MDispersedDT, as MassINI is used in, e.g., evaporation (which doesn’t vary by NBR_PARTIC). The results are shown in this notebook.

** usage of MassINI **
[3204] Me%Var%FMEvaporated = Me%Var%MEvaporated/Me%Var%MassINI
[3408] Me%Var%FMDispersed = (Me%Var%MDispersed+Me%Var%MBio)/Me%Var%MassINI
[4297] Me%Var%FMSedimented = Me%Var%MSedimented/Me%Var%MassINI
[4460] Me%Var%FMBio = Me%Var%MBio / Me%Var%MassINI
[4541] Me%Var%FMDissolved = Me%Var%MDissolved / Me%Var%MassINI

** Usage of MassOil **

[3169] Me%Var%MassOil = Me%Var%MassOil - (Me%Var%MEvaporatedDT) &
* Me%Var%DTOilInternalProcesses
[3363] Me%Var%MassOil = Me%Var%MassOil - (Me%Var%MdispersedDT) &
* Me%Var%DTOilInternalProcesses
[4279] Me%Var%MassOil = Me%Var%MassOil - (Me%Var%MSedimentedDT) &
* Me%Var%DTOilInternalProcesses

** MBio mass adjustment missing ** (I added it)

[4509] Me%Var%MassOil = Me%Var%MassOil - (Me%Var%MdissolvedDT) &
* Me%Var%DTOilInternalProcesses

MassINI is the Initial Mass of oil spill while MassOil is the mass of the oil that remains after being changed by weathering processes. I know this for sure b/c I wrote some print statements in ModuleOil_0D.F90. MassINI stays the same over time while MassOil decreases over time. The first set of print statements look like:

ModuleOil_0D.F90,  Me%Var%MassOil[3415]:    869632.212607649
 ModuleOil_0D.F90,  Me%Var%MassINI[3416]:    880175.032067522 

While the second hour of the run shows

ModuleOil_0D.F90,  Me%Var%MassOil[3415]:    761498.810614312
ModuleOil_0D.F90,  Me%Var%MassINI[3416]:    880175.032067522

etc.

So “FM” stands for Fraction of (initial) Mass.

[RachaelDMueller]

In this notebook, I show that the change in oil mass due to NBR_PARTIC is purely from the dispersion and biodegredation parameterizations. Changing NBR_PARTIC by a factor of 3 amplified dispersed mass by ~7-fold; whereas it changes overall mass by .14%. The ration between the Dispersed Mass / Total Mass = ~70 tonnes / 552 tonnes for the NBR_PARTIC = 1000 case. If Dispersed mass changes from 50 to 375 and total mass is 552 then we would see a greater change in the total mass from 2D thickness if MDispersed was being accounted for in the 2D thickness.

Need to take a careful look at MDispersed, Thickness_2D and the output to .sro file.

[SusanEAllen]

Susan's analysis shows that the results are moderately sensitive to the number of particles upto 10000 particles (20000 is very very similar to 10000) and that the run time is very similar upto 10000 particles with 20000 longer. So 10000 is the sweet spot.