Melt ponds on arctic sea-ice

globals[
  ;mean-ice-height is defined via the GUI
  ;smooth-cycles   is defined via the GUI
  density-ratio
  ;melt-rate-ice
  ;melt-rate-pond
  water-height-transition
  ;seepage-rate

  gravity
  kinematic-viscosity
  horizontal-permeability
  albedo-sea-water
  albedo-ice

  ;; physical parameters for the model
  time-step
  space-step ; not used for the moment

  ;utility variables
  max1
  max2
  initial-mean-ice
  time-to-half-of-initial-mean-ice
]

breed[drops drop]

; internal variables on the patches:
; "ice" is ice height in [cm]
; "water" is the melt water in [cm]
; "melt" represents the water recently melted in [cm]
; "albedo" the fraction of radiation reflected by the surface
patches-own[ice water albedo]

drops-own[water-content]

;; we affect the global physical parameters that are relevant for the model

to startup
  set time-step  0.5  ; expressed in days
  set space-step 100.  ; the lateral size of a patch expressed in cm

  ;set melt-rate-ice 1.2  ;in cm of ice per day
  ;set melt-rate-pond 2.0 ;in cm of ice per day, nominal value, check the literature
  set water-height-transition 10 ;in cm
  set density-ratio 0.8 ;ratio between ice and water mass densities , for less porous ice can be 0.9

  ;ifelse seepage? [
  ;set seepage-rate 0.8 ; in cm / day
  ;][
  ;set seepage-rate 0.0 ; in cm / day
  ;]

  set gravity 9.81 * 100 * (86400 * 86400 ) ;in  cm / day^2
  set kinematic-viscosity 1.e-6 * (100 * 100) * 86400  ;in cm^2 / day
  set horizontal-permeability  3.e-9 * (100 * 100); in cm^2

  set albedo-sea-water 0.1
  set albedo-ice 0.9

  set max1 0
  set max2 0
  set time-to-half-of-initial-mean-ice 0
end 


; generate a random smooth topography

to setup-topography
  if clear-previous-plots? [clear-all]
  startup
   ask patches
  [
  set ice random-float (2 * mean-ice-height)
  set water 0
  set albedo compute-albedo

  set pcolor scale-color cyan ice 0 (2 * mean-ice-height)
  ]
  ifelse(smooth-with-radius?) [
    smooth-ice-with-radius][
    smooth-ice-with-cycles]

  set initial-mean-ice compute-mean-ice

  reset-ticks
  tick ;this is to enable the plotting of the histogram of the initial configuration
  reset-ticks
end 

; this is the function to smooth the random ice field

to smooth-ice-with-cycles
  repeat smooth-cycles [
    ask patches [
      let sum-ice-neighbors sum [ice] of neighbors4
      let mean-ice (sum-ice-neighbors + ice) / 5
      set ice mean-ice + random-float 0.1
      color-field
    ]
  ]
end 

; this is the function to smooth the random ice field with a different procedure

to smooth-ice-with-radius
    ask patches [
      ;let sum-ice sum [ice] of patches in-radius smooth-radius
      ;let count-ice count patches in-radius smooth-radius
      ;let mean-ice (sum-ice  + ice) / ( count-ice + 1 )
      let mean-ice mean [ice] of patches in-radius smooth-radius
      set ice mean-ice
  ]
  ask patches[
      color-field
    ]
end 



; coloring functions

to color-field
  ifelse water > 0 [
    set pcolor scale-color blue (2. * mean-ice-height * density-ratio - water) 0 (2. * mean-ice-height * density-ratio)
  ][
  if ice > 0 [set pcolor scale-color grey ice -10 (2 * mean-ice-height)]
  ;same result as the following 2 lines
  ;let ice-color 80 + ((89.9 - 80) / (2 * mean-ice-height)) * ice
  ;set pcolor ice-color
    if ice = 0 [set pcolor turquoise] ;blue - 3]
  ]
end 



; melt ice

to melt-ice
  let actual-melted-volume 0
 if ice > 0 [
    ;; VERTICAL MELTING

    ;; 1) the following line implements conductive melting of ice
    if water = 0 [
      set actual-melted-volume  melt-rate-ice * time-step ;shall be corrected to include Stefan effect, see line below
      ;set actual-melted-volume  ( melt-rate-ice * time-step / ice * mean-ice-height )
    ]

    ;; 2) the following line implements the increased melt-rate for pondend ice : water enahnces melting (as in Luethje et al. paper)
    if water > 0 [
      ifelse water < water-height-transition[
      set actual-melted-volume (melt-rate-ice +  (melt-rate-pond - melt-rate-ice) * (water / water-height-transition)  )* time-step
      ][
      set actual-melted-volume melt-rate-pond * time-step
      ]
    ]

    ;; lateral melting
    if pond-lateral-melting?[set actual-melted-volume (actual-melted-volume + (melt-rate-pond * lateral-melting * time-step))]

    ;; Making melting happening
    if actual-melted-volume > ice [
      set actual-melted-volume ice  ;this is to avoid to melt more ice than what we have
      set water 0
    ]
    set ice (ice - actual-melted-volume ) ; melt has occurred
    sprout-drops 1 [
      set water-content (actual-melted-volume  * density-ratio)  ; the melted water is put into a pocket (a moving agent)
      ifelse pen-down? [pen-down][pen-up]
      hide-turtle
    ]
  ]
end 

;; seepage of meltwater

to seepage
  let seepage-amount (seepage-rate * time-step) ; the seepage amount per patch

    ask drops-here [
      ifelse water-content > seepage-amount[
        set water-content water-content - seepage-amount
        set seepage-amount 0
      ][
        set seepage-amount (seepage-amount - water-content)
        set water-content 0
        die
      ]
    ]


  if seepage-amount > 0 and water > 0[
    ifelse water > seepage-amount[
      set water (water - seepage-amount)
    ][
      set water 0
    ]
  ]
end 


;; move drops till conversion into water
;; flow assumes that the water produced by melt displaces instantaneously to position of minimum potential energy

to flow
  loop[
  let p min-one-of neighbors [ice + water]
   ifelse (ice + water) >  [ice + water] of p [
      move-to p
    ][
     set water (water + water-content) ; the pocket of melted water has reached a minimum height and it remains there
     if ice = 0 [set water 0]
     if not melt-ponds? [set water 0] ; remove all water if we are not interested in ponds
      die
  ]
  ]
end 

; move drops till conversion into water
;; flow2 assumes that the spped is finite. It depends on the local gradient of "ice + water" and by the size of the patch. Needs improvement.

to flow2
  let p min-one-of neighbors [ice + water]
  ;let displacement (gravity / kinematic-viscosity) * (((ice + water) - [ice + water] of p) / space-step) * time-step
  ;let gradient (((ice + water) - [ice + water] of p) / space-step)
  let h atan space-step ((ice + water) - [ice + water] of p)
  let angle (90 - h) mod 360
  let horizontal-velocity  2 ;1000 * space-step  * sqrt (abs (gravity * (sin angle))  / (2 * space-step) )
  let displacement (horizontal-velocity * time-step)
   ifelse displacement >  distance p [
      move-to p
    ][
     set water (water + water-content) ; the pocket of melted water has reached a minimum height and it remains there
     if ice = 0 [set water 0]
     if not melt-ponds? [set water 0] ; remove all water if we are not interested in ponds
     die
  ]
end 

to refreeze
  ask drops[
    set water (water + water-content) ; the pocket of melted water has reached a minimum height and it remains there
    if ice = 0 [set water 0]
    if not melt-ponds? [set water 0] ; remove all water if we are not interested in ponds
    die
  ]
  ask patches[
  set ice ice + (water / density-ratio)
  set water 0
  color-field
  ]
  tick
end 

; a function for patches to account for the contribution to melting due to the presence of nearby ponds

to-report lateral-melting
  let ice-here ice
  let melt-volume 0
  ask neighbors4 with [ice < ice-here][
  ;ask neighbors with [ice < ice-here][
    ifelse (ice + water > ice-here)[
      set melt-volume (melt-volume + (ice-here - ice))
    ][
      set melt-volume (melt-volume + water)
    ]
  ]
  set melt-volume (melt-volume / space-step)
  report melt-volume
end 


;; MAIN
;; this is the procedure for the loop over time

to melt-and-flow
  if pen-down? [cd] ; to clear previous drawings

  ask patches[
    ;; melt the ice
    melt-ice
    ;; seepage
    if seepage? [seepage]
  ]

  ask drops[
    ;; move water
    if water-flowing-mode = 1 [flow]
    if water-flowing-mode = 2 [flow2]
  ]

  ask patches[
    ;; recoloring the map
    color-field
    ;; estimate albedo
    set albedo compute-albedo
  ]

  if (count patches with [ice = 0] = count patches )[
    stop
  ]
  tick
end 

; reports used in plots

to-report compute-mean-ice
  report  mean [ice] of patches
end 

to-report compute-mean-water
  report  mean [water] of patches
end 

to-report compute-max-mean-water
  let water-depth mean [water] of patches
  if water-depth > max2 [ set max2 water-depth]
  report max2
end 

to-report compute-std-ice
  let avg compute-mean-ice
  let var (mean [ice * ice] of patches)
  report sqrt (var - ( avg * avg ))
end 

; we need to improve this, it's just for test

to-report compute-albedo
  let alpha albedo-sea-water  ; we start from total absorption (like if there is sea everywhere)
  if ice > 0 [
  ifelse water > 0 [
    let kappa (-1 / water-height-transition) * ln ( albedo-sea-water / albedo-ice )
    set alpha albedo-ice * exp (- water / water-height-transition) ; albedo coeff of pond
    ][
    set alpha albedo-ice ; albedo coeff of ice
    ]
  ]
  report alpha
end 

to-report max-pond-coverage
  let value 100 * (count patches with [water > 0]) / count patches
  if value > max1 [set max1 value]
  report max1
end 

to-report time-to-50%-ice-melt
  if compute-mean-ice <  0.5 * initial-mean-ice and time-to-half-of-initial-mean-ice = 0 [
    set time-to-half-of-initial-mean-ice (ticks * time-step)
  ]
    report time-to-half-of-initial-mean-ice
end