Hyporheic Exchange Model: Interface

breed [ particles particle ]

particles-own [ residence-time in-bank? ]
globals [Boundary constant maxflux k2 height time-counter time-dx HydrCond2 hm2 Boundary2 xcutoff ChannelVelocity2 depth2]

to-report surface-function [ x ]
  report height * sin ( k2 * x ) + Slope * x
end 

to setup
  clear-all
 set constant 240 / (3 * lambda)

 set k2 ( 360 / (constant * lambda ))
 set Boundary BedDepth
 set height (( BedformHeight / 2 ) * constant )
 set time-counter 0
 set time-dx 100
 set HydrCond2 constant * HydrCond
 set ChannelVelocity2 ChannelVelocity * constant
 set depth2 depth * constant
 ifelse (( height * 2 / depth2 ) > .34)
  [set hm2 .28 * ( ( ChannelVelocity2 ) ^ 2 ) / ( 2 * 9.8 * constant * 100  ) * (( height * 2 / depth2 ) / .34 ) ^ ( 3 / 2 ) ]
  [set hm2 .28 * ( ( ChannelVelocity2 ) ^ 2 ) / ( 2 * 9.8 * constant * 100  ) * (( height * 2 / depth2 ) / .34 ) ^ ( 3 / 8 ) ]

 set Boundary2 constant * Boundary
 set xcutoff constant * lambda * floor( 240 / (constant * lambda))

 set maxflux 1.5 * ( max [flux-weighting pxcor] of patches with [pycor = 0] )
 ;show maxflux

 ask patches
 [
  ifelse ( pycor < surface-function (pxcor))
    [ set pcolor brown + 2 ] ; sand color
    [ set pcolor blue - 2 ]  ; water color

 if ( pycor < -1 * boundary2 )
    [ set pcolor black ]

 if ( pxcor > xcutoff )
    [ set pcolor black ]


   if ( pxcor < 0 )
    [ set pcolor black ]
  ]

 repeat NumParticles
 [
  drop-particle
 ]
do-plot
reset-ticks
end 

to-report flux-weighting [ x ]
let exp1 1
let exp2 exp ( k2 * (2 * pi / 360 ) * (-2 * boundary2 ) )
let expsum exp1 + exp2

let yvel  ((-1 * (cos ( x * k2 )) * ( ( exp1 - exp2) * (k2 * 2 * pi / 360 )) * hm2)) * HydrCond2 / porosity ;
let xvel  ((k2 * (2 * pi / 360 ) * ( sin ( x * k2 )) * ( expsum )  * hm2 ) ) * HydrCond2 / porosity;
let derivative (( (1) * height * 2 * pi / (constant * lambda )) * cos ( ( x) * k2 ) )
report derivative * xvel - yvel
end 

to-report random-flux-weighted-x

  let x random-float xcutoff
  let threshold random-float maxflux   ; number bigger than any flux-weighting value
  while [ threshold > flux-weighting x ]
  [
     set x random-float xcutoff
     set threshold random-float maxflux
  ]
  report x
end 

to drop-particle
 create-particles 1 [
   set xcor random-flux-weighted-x
   set ycor surface-function xcor
   set in-bank? true
   set residence-time 0
   set size .5
   set color green
   set shape "circle"
   stamp
   set size 2

   if (show-paths?) [ pd ]
 ]
end 

to go
 ask particles with [ in-bank? ]
 [
  move
 ]

 if re-drop?
 [
 while [ count particles with [ in-bank? ] < NumParticles ]
 [
  drop-particle
 ]
]
 set time-counter time-counter + time-dx
if (not re-drop?)
 [
  do-plot
 ]
tick
end 

to move

let ycor2 ycor - surface-function( xcor )

let exp1 exp ( k2 * (2 * pi / 360 ) * ycor2 )
let exp2 exp ( k2 * (2 * pi / 360 ) * (-2 * boundary2 - ycor2 ) )
let exp3 exp ( k2 * (2 * pi / 360 ) * (-2 * boundary2) )
let expsum exp1 + exp2

let yvel  ((-1 * (cos ( xcor * k2 )) * ( ( exp1 - exp2) * (k2 * 2 * pi / 360 )) * hm2)) / (1 + exp3) * HydrCond2 / porosity ;
let xvel  ((k2 * (2 * pi / 360 ) * ( sin ( xcor * k2 )) * ( expsum )  * hm2 ) ) / (1 + exp3) * HydrCond2 / porosity;

 set xcor xcor +  time-dx * (xvel - HydrCond2 * Slope / (porosity * (1 + Slope ^ 2))) ; / constant

 if (ycor + yvel > max-pycor)
 [ show "big problem -- particle lost above screen" die ]
 if (ycor + yvel < min-pycor)
 [
  ;show "Oh no, I'm lost!"
  die
 ]
  set ycor ycor + time-dx * (yvel - (HydrCond2 * Slope ^ 2) / (porosity * (1 + Slope ^ 2))); + HydrCond2 / porosity );/ constant

 set residence-time residence-time + time-dx

if xcor >= xcutoff
  [ pen-up

   set ycor ycor - Slope * xcutoff
    set xcor xcor - xcutoff

  if show-paths?
    [pen-down]

  ]

  if xcor <= 0
  [ pen-up

      set ycor ycor + Slope * xcutoff

    set xcor xcor + xcutoff

  if show-paths?
    [pen-down]
  ]

 set in-bank? (ycor <= surface-function xcor)
 ifelse (in-bank?)
 [
  set color blue
 ]
 [
  set color magenta
  ; show residence-time
  set size 0.5
  ; we want them to stamp them right on the function, for better visuals
  set ycor surface-function xcor
  stamp
  ht
 ]
end 

to do-plot
 if (count particles != 0)
[
  set-current-plot "Cumulative Residence Time"

if time-counter > time-dx
  [set-plot-x-range 0 1 + round (log (1 + Median [residence-time] of particles) 10 ) * 3]
  ; keep increasing y-range, proportionally to number of particles
  ;set-plot-y-range 0 (ceiling (count particles / 1900))
  set-plot-y-range 0 1
  ;set-histogram-num-bars 54
  ;histogram [ residence-time ] of particles with [ not in-bank? ]
  plotxy (log (( time-counter  + 1 )) 10)  ( count particles with [ in-bank? ] / count particles )

]
end 

;; COMMENTED OUT 10/15/2021 to allow running in NetLogoWeb
;to export-residence-times
; let file user-new-file
; if ( file != false )
; [
;  if ( file-exists? file )
;    [ file-delete file ]
;  file-open file
; ]
; let res-time-list [ residence-time ] of particles with [ not in-bank? ]
; foreach res-time-list
; [ ?1 ->
;  file-type ?1
;  file-type "\r\n"
; ]
; file-close
;
;end