; ================================================================================
; SINKHOLE SIMULATION - Educational Version
; ================================================================================
; Version: v7.0
; Description: Urban area sinkhole formation mechanism simulation
; 
; Author: Park, Hongjoon¹
; Affiliation: 
;   ¹ Chonnam National University, Adjunct Professor
; 
; Copyright (c) 2025 Park, Hongjoon. All rights reserved.
; 
; License: CC BY-NC-SA 4.0
;   (Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International)
;   
;   You are free to:
;   - Share and adapt this work for educational and research purposes
;   - Must give appropriate credit and indicate if changes were made
;   - NonCommercial use only
;   - Derivative works must use the same license
;   
;   For commercial use or other permissions, contact: [your.email@jnu.ac.kr]
;
; Citation: If you use this model in your research, please cite:
;   Park, H. (2025). Sinkhole Simulation Model v7.0. 
;   Chonnam National University.
;   [Will be updated with publication information]
;
; Development Period: 2025-2026
; Last Updated: January 2026
; ================================================================================

; ================================================================================
; 1. Global Variables and Constants Definition
; ================================================================================

globals [ 
; Statistical variables
total-soil-loss          ; Total soil loss amount
cavity-size              ; Total cavity size
collapse-occurred?       ; Whether collapse occurred
time-to-collapse         ; Time taken until collapse
groundwater-flow         ; Groundwater flow velocity
effective-stress         ; Effective stress
danger-index             ; Danger index (0-100)
excavation-position      ; Construction site position
show-cavities?           ; Show/hide cavities

; Threshold constants
WATER-FLOW-THRESHOLD        ; Water flow generation threshold
EXCAVATION-EROSION-THRESHOLD ; Excavation erosion threshold
]

breed [soil-particles particle]
breed [water-flows flow]
breed [cavities cavity]
breed [excavations excavation]
breed [persons person]

soil-particles-own [
soil-type                ; "sand", "clay", "gravel"
is-fine-particle?        ; Whether fine-grained soil
stability                ; Stability (0-100)
lost?                    ; Whether lost
]

water-flows-own [
velocity                 ; Flow velocity
carrying-capacity        ; Carrying capacity
]

cavities-own [
cavity-depth             ; Cavity depth
cavity-width             ; Cavity width
growth-rate              ; Growth rate
]

patches-own [
layer-type               ; Layer type
has-pipe?                ; Pipe existence
pipe-health              ; Pipe health (0-100, higher is better)
water-saturation         ; Groundwater saturation (0-100)
soil-density             ; Soil density
permeability             ; Permeability (0-100)
has-cavity?              ; Cavity existence
cavity-volume            ; Cavity volume
]

; ================================================================================
; 2. Simulation Initialization
; ================================================================================

to setup
clear-all

; ============================================================
; Define simulation constants (modify only here when needed)
; ============================================================
set WATER-FLOW-THRESHOLD 30           ; Minimum saturation for water flow generation (%)
set EXCAVATION-EROSION-THRESHOLD 35   ; Minimum saturation for excavation erosion (%)
; ============================================================

; Initialize global variables
set total-soil-loss 0
set cavity-size 0
set collapse-occurred? false
set time-to-collapse 0
set danger-index 0
set excavation-position random 15 - 10
set show-cavities? true

; Setup soil layer structure
setup-soil-layers

; Setup underground pipes
setup-underground-pipes

; Setup groundwater
setup-groundwater

; Setup soil particles
setup-soil-particles

reset-ticks
end 

; Procedure for cavity display toggle

to cavity-display
  set show-cavities? not show-cavities?
  ifelse show-cavities? [
    ask cavities [ show-turtle ]
  ][
    ask cavities [ hide-turtle ]
  ]
end 

; ================================================================================
; 3. Soil Layer Structure Setup
; ================================================================================

to setup-soil-layers
; Define layer boundaries
let sky-bottom (max-pycor - 9)
let grass-layer (max-pycor - 9)
let surface-bottom (max-pycor - 16)
let fills-bottom (max-pycor - 22)
let weathered-soil-bottom (max-pycor - 30)
let weathered-rock-bottom (max-pycor - 38)

create-persons surface-load [
  setxy random-xcor grass-layer + 1
  set size 1.5 
  set shape "person"
  ]

ask patches [
 ; Determine layer type and color based on y-coordinate of each patch
 if pycor > sky-bottom [
   set layer-type "sky"
   set pcolor sky
 ]

 if pycor = grass-layer [
   set layer-type "grass"
   set pcolor green
 ]

 if pycor <= grass-layer - 1 and pycor > surface-bottom [
   set layer-type "surface"
   set pcolor brown + 2
 ]

 if pycor <= surface-bottom and pycor > fills-bottom [
   set layer-type "fills"
   set pcolor brown
 ]

 if pycor <= fills-bottom and pycor > weathered-soil-bottom [
   set layer-type "weathered-soil"
   set pcolor grey + 2
 ]

 if pycor <= weathered-soil-bottom and pycor > weathered-rock-bottom [
   set layer-type "weathered-rock"
   set pcolor grey
 ]

 if pycor <= weathered-rock-bottom [
   set layer-type "rock-layer"
   set pcolor grey - 2
 ]

 ; Set basic soil properties
 set soil-density compaction-degree
 set permeability 100 - compaction-degree
 set has-cavity? false
 set cavity-volume 0
 set has-pipe? false
]
create-excavations 1 [
  setxy excavation-position 45 - 7
  set size 3
  set shape "flag"
  set color grey + 2
]
end 

; ================================================================================
; 4. Underground Pipe Setup (Cluster-based)
; ================================================================================

to setup-underground-pipes
; Cluster-based pipe placement system

; Generate 2-4 cluster centers
let num-clusters 2 + random 3
let cluster-centers []  ; Create a list variable

repeat num-clusters [
 ; Select random position within fill layer as cluster center
 let center-patch one-of patches with [layer-type = "fills"]
 if center-patch != nobody [
   set cluster-centers lput center-patch cluster-centers
 ]
]

; Place pipes densely around each cluster center
foreach cluster-centers [ center ->
 ask patches with [layer-type = "fills"] [
   ; Calculate distance from cluster center
   let dist distance center

   ; Place pipes with probability inversely proportional to distance (higher probability closer)
   let placement-prob 0

   if dist < 5 [
     ; Very close distance: 80% probability
     set placement-prob 80
   ]
   if dist >= 5 and dist < 10 [
     ; Medium distance: 50% probability
     set placement-prob 50
   ]
   if dist >= 10 and dist < 15 [
     ; Far distance: 20% probability
     set placement-prob 20
   ]

   ; Place pipe according to probability
   if random 100 < placement-prob [
     set has-pipe? true
     set pipe-health pipe-condition  ; Changed: get value from slider
     set pcolor black
   ]
 ]
]
end 

; ================================================================================
; 5. Groundwater Setup
; ================================================================================

to setup-groundwater
let grass-depth 10
let surface-depth 6
; Set initial groundwater depth from below pipe layer (6~29)
let groundwater-base-depth (grass-depth + surface-depth + initial-groundwater-depth)

ask patches [
 ifelse pycor < max-pycor - groundwater-base-depth [
   set water-saturation random 50  ; Set groundwater saturation randomly below 50
   set pcolor pcolor + (random-float 1.5)
 ] [
   set water-saturation 0 ; Areas outside groundwater depth have 0 saturation
 ]
]
end 

; ================================================================================
; 6. Soil Particle Setup
; ================================================================================

to setup-soil-particles
ask patches with [ layer-type = "surface"
                  or layer-type = "weathered-soil"
                  or (layer-type = "fills" and has-pipe? = false) ][
 ; Soil particle generation probability: 80%
 if random 100 < 80 [
   sprout-soil-particles 4 [
     set shape "circle"
     set size 0.6
     set lost? false
  
     ; Set properties by soil type
     if soil-type-selection = "sand" [
       set soil-type "sand"
       set is-fine-particle? true
       set color yellow
       set stability 20 + random 30
     ]
  
     if soil-type-selection = "gravel" [
       set soil-type "gravel"
       set is-fine-particle? false
       set color violet
       set stability 40 + random 40
     ]
  
     if soil-type-selection = "clay" [
       set soil-type "clay"
       set is-fine-particle? true
       set color orange
       set stability 60 + random 30
     ]
  
     if soil-type-selection = "mixed" [
       let rand random 100
       ifelse rand < 33 [
         set soil-type "sand"
         set is-fine-particle? true
         set color yellow
         set stability 20 + random 30
       ] [
         ifelse rand < 66 [
           set soil-type "clay"
           set is-fine-particle? true
           set color orange
           set stability 60 + random 30
         ] [
           set soil-type "gravel"
           set is-fine-particle? false
           set color violet
           set stability 40 + random 40
         ]
       ]
     ]
  
     ; Adjust stability according to compaction degree
     set stability stability * (compaction-degree / 100)
   ]
 ]
]
end 

; ================================================================================
; 7. Main Simulation Loop
; ================================================================================

to go
if collapse-occurred? [ 
  ask patch 0 43 [ set plabel "Finished"]
  stop 
]

; Rainfall effect
apply-rainfall

; Update groundwater
update-groundwater

; Create water flow
create-water-flow

; Pipe leakage erosion and cavity formation
pipe-erosion-and-cavity

; Excavation effect
excavation-effect

; Fine particle erosion
erode-particles

; Cavity expansion
expand-cavities

; Surface load
apply-load

; Check collapse
check-collapse

; Update statistics
update-statistics

tick
if not collapse-occurred? [
 set time-to-collapse time-to-collapse + 1
]
end 

; ================================================================================
; 8. Rainfall Effect
; ================================================================================

to apply-rainfall
; rainfall-intensity: limited to 0~50 by slider
; water-saturation: 0~100
if rainfall-intensity > 0 [ ; Execute only when above 0
 ; Sky visual effect
 ask patches with [layer-type = "sky"] [
   if random 100 < rainfall-intensity [
     set pcolor sky - random-float 0.5
     set pcolor sky + random-float 0.8
   ]
 ]

 ; Increase groundwater saturation
 ask patches with [water-saturation > 0] [
   set water-saturation min list 100 (water-saturation + rainfall-intensity / 50)
 ]

 ; Additional groundwater saturation from rainwater infiltration, up to 4 with rainfall intensity probability
 ask patches with [layer-type = "surface" or layer-type = "fills"] [
   if water-saturation < 90 and random 100 < rainfall-intensity [
     set water-saturation water-saturation + random 5
   ]
 ]
]

; Clear sky when no rainfall
if rainfall-intensity = 0 [
 ask patches with [layer-type = "sky"] [
   set pcolor sky
 ]
]
end 

; ================================================================================
; 9. Groundwater Update: water-saturation changes
; ================================================================================

to update-groundwater
; Groundwater diffusion: increase water-saturation values in 4 neighbors
; groundwater-change-rate: 0~30
ask patches with [water-saturation > 0] [ ; When groundwater saturation value is not 0
 ask neighbors4 with [water-saturation < 100] [
   set water-saturation min list 100 (water-saturation + groundwater-change-rate / 30)
 ]
]

; Groundwater visualization
; Excluding sky and grass layers, color change range (0.3), applied (30% probability)
ask patches with [water-saturation > 50 and
                 layer-type != "sky" and
                 layer-type != "grass"] [
 if random 100 < 30 [
   set pcolor pcolor + (random-float 0.6 - 0.3)  ; Range -0.3 to +0.3
 ]
]

; Drainage due to excavation, typical excavation depth is 0~25m
; excavation-depth: 0~25
if excavation-depth > 0 [
 let drainage-zone-top 36 - excavation-depth ; Subtracting 0 from 36 starts at surface layer
 ; Water drains from patches above excavation depth with saturation over 50
 ask patches with [pycor > drainage-zone-top and water-saturation > 50] [
   set water-saturation max list 0 (water-saturation - random-float 5.0)
 ]
]
end 

; ================================================================================
; 10. Water Flow Creation
; ================================================================================

to create-water-flow
ask water-flows [ die ]

; Using constant: WATER-FLOW-THRESHOLD ← currently 30
; groundwater-change-rate: 0~30
; permeability: initial value is 100 - compaction-degree in layer structure setup
ask patches with [water-saturation > WATER-FLOW-THRESHOLD and has-cavity?] [
 if random 100 < 50 [
   sprout-water-flows 1 [
     set shape "arrow"
     set color cyan
     set size 1.0
     set velocity groundwater-change-rate * permeability / 100
     set carrying-capacity velocity * 2
  
     ; Set flow direction: cavities have highest priority, pipes with low health are next
     if any? patches with [has-pipe? and pipe-health < 50] [
       face one-of patches with [has-pipe? and pipe-health < 50]
     ]
     if any? patches with [has-cavity?] [
       face one-of patches with [has-cavity?]
     ]
   ]
 ]
]

; Water flow movement
ask water-flows [
 forward velocity / 10
 if not can-move? 1 [ die ]
]
end 

; ================================================================================
; 11. Pipe Leakage Erosion and Cavity Formation
; ================================================================================

to pipe-erosion-and-cavity
  ; pipe-health is initialized with pipe-condition
  
  ask patches with [has-pipe?] [
    let erosion-prob 0
    
    ; 1) Base erosion rate (reflecting pipe-health): probability determined by multiplying 100-pipe-health by coefficient (0.15)
    let base-erosion (100 - pipe-health) * 0.15
    
    ; 2) Weight calculation (convert each variable's contribution to ratio between 0-1)
    ; Increases up to 1.5x as water-saturation increases
    let water-factor 1 + (water-saturation / 100 * 0.5)
    ; Increases up to 1.3x as rainfall-intensity increases
    let rain-factor 1 + (rainfall-intensity / 50 * 0.3)
    ; Increases up to 1.2x as excavation-depth increases
    let excavation-factor 1 + (excavation-depth / 25 * 0.2)
    
    ; 3) Determine final probability
    set erosion-prob base-erosion * water-factor * rain-factor * excavation-factor
    
    ; Erosion occurrence --> Fine particle erosion: lost? = true
    ; Non-linear increase --> P/100 * P/100 = P^2 / 10000 (probability overlaps twice)
    ; Loss speed increases exponentially even with slight increase in erosion-prob
    if random 100 < erosion-prob [      
      ask soil-particles-here [
        if is-fine-particle? and random 100 < erosion-prob [
          set lost? true
          set total-soil-loss total-soil-loss + 1
          die
        ]
      ]
      ; Set lower erosion probability for surrounding area (neighbors) (* 0.4)
      ask neighbors [
        if random 100 < (erosion-prob * 0.4) [
          ask soil-particles-here [
            if is-fine-particle? and random 100 < (erosion-prob * 0.4) [
              set lost? true
              set total-soil-loss total-soil-loss + 1
              die
            ]     
          ]
        ]
      ]
        
    ; Cavity formation: if soil particles = 0
    ; → has-cavity? true, cavity-volume 0.5, create cavity
    
    if not has-cavity? and count soil-particles-here = 0 [
      set has-cavity? true
      set cavity-volume 0.5
      set pcolor black
      
      sprout-cavities 1 [
        set shape "circle"
        ;set color red
        set color [255 0 0 128]
        set size 1
        set cavity-depth 1
        set cavity-width 1
        set growth-rate 0.1
      ]
    ]
  ]    
]
end 

; ================================================================================
; 12. Excavation Effect (soil particle outflow, cavity formation)
; ================================================================================

to excavation-effect
let excavation-zone-top max-pycor - 10 ; 35
let excavation-zone-bottom excavation-zone-top + 1 - excavation-depth  ; When 1, equals 35

; Excavation position (excavation-position)
if excavation-depth > 5 [  ; Generate only when excavation depth exceeds 5
 ask patches with [pycor <= excavation-zone-top and
                   pycor >= excavation-zone-bottom and
                   pxcor = excavation-position and
                   water-saturation > EXCAVATION-EROSION-THRESHOLD][ 
    let excavation-erosion-prob (excavation-depth / 25 * 0.01) ; Loss occurs with low probability
    ask soil-particles-here [
     ; Soil particle loss due to excavation
     if random-float 1.0 < excavation-erosion-prob [
       set lost? true
       set total-soil-loss total-soil-loss + 1
       die
     ]
    ]
   ; Cavity formation due to excavation
   if not has-cavity? and count soil-particles-here = 0 [
     set has-cavity? true
     set cavity-volume 0.5
     set pcolor black
  
     sprout-cavities 1 [
       set shape "circle"
       set  color [255 0 0 128]
       set size 1
       set cavity-depth 1
       set cavity-width 1
       set growth-rate 0.1
     ]
   ]
 ]
]
end 

; ================================================================================
; 13. Fine Particle Erosion (groundwater flow process -> fine particle erosion)
; ================================================================================

to erode-particles
ask water-flows [
 ask soil-particles in-radius 2 [
   if is-fine-particle? [
     ; Calculate erosion probability
     let base-erosion-prob ([velocity] of myself * [permeability] of patch-here / 100)
     let erosion-prob base-erosion-prob * (100 - stability) / 100
  
     if random-float 1.0 < erosion-prob / 50 [
       set lost? true
       set total-soil-loss total-soil-loss + 1
       die
     ]
   ]
 ]
]
end 

; ================================================================================
; 14. Cavity Expansion
; ================================================================================

to expand-cavities
ask cavities [
 ; Growth rate: 0.05 + groundwater change rate/300 + total soil loss/2000
 set growth-rate 0.05 + (groundwater-change-rate / 30 / 100) + (total-soil-loss / 2000)

 ; Increase cavity size
 set cavity-depth cavity-depth + growth-rate
 set cavity-width cavity-width + growth-rate / 3 ;2
 set size cavity-width

 ; Update cavity volume
 ask patch-here [
   set cavity-volume cavity-volume + [growth-rate] of myself

   ; Horizontal expansion: cavity-volume > 15, probability less than 15
   if cavity-volume > 15 [
     ask neighbors4 [
       if not has-cavity? and random 100 < 15 [
         set has-cavity? true
         set cavity-volume 0.5
         set pcolor black
      
         sprout-cavities 1 [
           set shape "circle"
           set color [255 0 0 128]
           set size 0.8
           set cavity-depth 1
           set cavity-width 1
           set growth-rate 0.1
         ]
       ]
     ]
   ]
 ]

 ; Vertical expansion: probability less than 25
 if random 100 < 25 [
   let target-patch patch-at 0 1

   if target-patch != nobody [
     ask target-patch [
       if not has-cavity? [  ; Soil particle erosion (probability less than 35)
         ask soil-particles-here [
           if random 100 < 35 [
             set lost? true
             set total-soil-loss total-soil-loss + 1
             die
           ]
         ]
      
         ; Cavity formation (< 2), volume reduction (/3)
         if count soil-particles-here < 2 [
           set has-cavity? true
           set cavity-volume [cavity-volume] of myself / 3
           set pcolor black
         ]
       ]
     ]
   ]
 ]
]

; Calculate total cavity size
set cavity-size sum [cavity-volume] of patches with [has-cavity?]
end 

; ================================================================================
; 15. Surface Load Application
; ================================================================================

to apply-load
ask patches with [layer-type = "surface"] [
 ; When there is a large cavity within radius 3
 if any? patches in-radius 3 with [has-cavity? and cavity-volume > 3] [
   if random 100 < (surface-load / 2) [
     let below-patch patch-at 0 -1
  
     if below-patch != nobody [
       ask below-patch [
         ask soil-particles-here [
           set lost? true
           set total-soil-loss total-soil-loss + 1
           die
         ]
       ]
     ]
   ]
 ]
]
end 

; ================================================================================
; 16. Collapse Check
; ================================================================================

to check-collapse
ask patches with [layer-type = "surface"] [
 let cavity-below patches in-radius 4 with [has-cavity? and cavity-volume > 15]
 if any? cavity-below [
   ; Calculate support strength: coefficient → /5
   let support-strength soil-density * (count soil-particles in-radius 3) / 5

   ; Calculate load pressure: coefficient → /20
   let load-pressure surface-load + cavity-size / 20

   ; Collapse occurs when support strength < load
   if support-strength < load-pressure [
     display  ; Force display of screen changes so far
     wait 0.05 ; Pause for 0.05 seconds

     ; Find position (x, y) of largest cavity
     let largest-cavity max-one-of cavity-below [cavity-volume]
     let collapse-x [pxcor] of largest-cavity
     let collapse-y [pycor] of largest-cavity
  
     ; Calculate dynamic size proportional to cavity size
     let total-cavity-volume sum [cavity-volume] of cavity-below
     let size-factor total-cavity-volume / 20
  
     ; Minimum/maximum size limits
     let collapse-width max list 3 (min list 8 (3 + size-factor * 0.5))
     let collapse-depth max list 4 (min list 10 (4 + size-factor * 0.6))
  
     ; Determine collapse type
     ifelse compaction-degree < 50 or groundwater-change-rate > 15 [
       ; Punching collapse (narrow and deep) - based on cavity position
       ask patches with [pycor > collapse-y - collapse-depth and
                        pycor < max-pycor - 8 and
                        abs(pxcor - collapse-x) < collapse-width] [
         if layer-type != "sky" [
           set pcolor black
           set has-cavity? true
           ask soil-particles-here [ die ]
         ]
       ]
    
       ; Surface subsidence area
       ask patch collapse-x (max-pycor - 10) [
         ask patches in-radius (collapse-width * 0.7) [
           if layer-type != "sky" [
             set pcolor red - 2
           ]
         ]
       ]
     ] [
       ; Trough collapse (wide and shallow) - based on cavity position
       ask patches with [pycor > collapse-y - (collapse-depth * 0.7) and
                        pycor < max-pycor - 8 and
                        abs(pxcor - collapse-x) < (collapse-width * 1.3)] [
         if layer-type != "sky" and random 100 < 60 [
           set pcolor black
           set has-cavity? true
           ask soil-particles-here [
             if random 100 < 70 [ die ]
           ]
         ]
       ]
    
       ; Surface subsidence area
       ask patch collapse-x (max-pycor - 10) [
         ask patches in-radius (collapse-width * 1.0) [
           if layer-type != "sky" [
             set pcolor red + 2
           ]
         ]
       ]
     ]
      set collapse-occurred? true
   ]
 ]
]
end 

; ================================================================================
; 17. Statistics Update
; ================================================================================

to update-statistics
; Calculate danger index
let cavity-factor min list 50 (cavity-size * 2)
let soil-loss-factor min list 30 (total-soil-loss / 10)
let load-factor min list 20 (surface-load / 5)

set danger-index cavity-factor + soil-loss-factor + load-factor

if collapse-occurred? [
 set danger-index 100
]

; Average groundwater flow velocity
ifelse any? water-flows [
 set groundwater-flow mean [velocity] of water-flows
] [
 set groundwater-flow 0
]

; Effective stress
set effective-stress compaction-degree *
                    (1 - cavity-size / 100) *
                    (1 - total-soil-loss / 1000)
end 

; ================================================================================
; End
; ================================================================================