tReorder figures - cngf-pf-exp1 - experiments for first paper with continuum granular model
(HTM) git clone git://src.adamsgaard.dk/manus_continuum_granular1_exp
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) Submodules
(DIR) README
(DIR) LICENSE
---
(DIR) commit ab7a19a2bc822b9130649f6a99492b68e2fd7bf1
(DIR) parent 49bd09b37c66b776b58228e0154336f5efab996b
(HTM) Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Mon, 8 Jul 2019 17:50:01 +0200
Reorder figures
Diffstat:
M fig2/Makefile | 61 ++++++++++++++++---------------
M fig2/fig.gp | 47 +++++++++++--------------------
M fig3/Makefile | 5 +----
M fig3/fig.gp | 54 ++++++++++++++-----------------
M fig4/Makefile | 35 +++++++++++++++++++++++++++----
M fig4/fig.gp | 84 ++++++++++++-------------------
M fig5/Makefile | 3 +--
M fig5/fig.gp | 11 +++++------
M fig6/Makefile | 14 +++++++++++++-
M fig6/fig.gp | 95 ++++++++++++++++++-------------
M fig7/Makefile | 38 +------------------------------
M fig7/fig.gp | 55 ++++++++++++++++++++++---------
12 files changed, 251 insertions(+), 251 deletions(-)
---
(DIR) diff --git a/fig2/Makefile b/fig2/Makefile
t@@ -1,42 +1,43 @@
BIN = ../1d_fd_simple_shear/1d_fd_simple_shear
FIG = $(shell basename $(CURDIR))
+RESOLUTION := 1000
+LENGTH := 1.0
+NORMAL_STRESS := 85e3
+
default: ../$(FIG).pdf
-sim.output00000.txt: $(BIN)
+two_rivers.txt: $(BIN) Makefile
/bin/sh -c '\
- ./$< --resolution 200 --length 8.0 --normal-stress 150e3 \
- --fluid --fluid-permeability 2e-17 --fluid-pressure-top 50e3 \
- --fluid-pressure-ampl 50e3 \
- --fluid-pressure-freq $$( echo "1.0/(3600*24)" | bc -l ) \
- --file-interval $$( echo "60*10" | bc -l ) \
- --time-end $$( echo "3600*24*7" | bc -l ) sim'
-
-timeseries.txt: sim.output00000.txt
+ ./$< --resolution $(RESOLUTION) \
+ --length $(LENGTH) \
+ --normal-stress $(NORMAL_STRESS) \
+ --friction-coefficient 0.321 --cohesion 14.0e3 \
+ --stress-ratio 0.5 \
+ --grain-size 0.1e-3 \
+ --normalize' > $@
+
+storglaciaren.txt: $(BIN) Makefile
/bin/sh -c '\
- for f in sim.output*.txt; do \
- tail -n 1 "$$f" | cut -f2- >> $@; \
- done'
+ ./$< --resolution $(RESOLUTION) \
+ --length $(LENGTH) \
+ --normal-stress $(NORMAL_STRESS) \
+ --friction-coefficient 0.494 --cohesion 5.0e3 \
+ --stress-ratio 0.594 \
+ --grain-size 0.1e-3 \
+ --normalize' > $@
-max_strainrate_depth.txt: sim.output00000.txt
+upb.txt: $(BIN) Makefile
/bin/sh -c '\
- for f in sim.output*.txt; do \
- awk " \
- BEGIN{getline; depth=8.0; max=0.0} \
- NF{ \
- if (\$$6 >= max && \$$6+0 > 1e-300) \
- { \
- max=\$$6; \
- depth=\$$1; \
- } \
- } \
- END{ \
- print depth,\"\t\",max \
- }" \
- "$$f" >> $@; \
- done'
-
-../$(FIG).pdf: fig.gp timeseries.txt max_strainrate_depth.txt
+ ./$< --resolution $(RESOLUTION) \
+ --length $(LENGTH) \
+ --normal-stress $(NORMAL_STRESS) \
+ --friction-coefficient 0.443 --cohesion 3.0e3 \
+ --stress-ratio 0.543 \
+ --grain-size 0.1e-3 \
+ --normalize' > $@
+
+../$(FIG).pdf: fig.gp two_rivers.txt storglaciaren.txt upb.txt
gnuplot $< > $@
clean:
(DIR) diff --git a/fig2/fig.gp b/fig2/fig.gp
t@@ -1,36 +1,23 @@
#!/usr/bin/env gnuplot
-reset
-set terminal pdfcairo color size 7.5 cm, 8.5 cm
+reset
-set multiplot layout 3,1 \
- margins 0.19,0.85,0.23,0.95 \
- spacing 0.03,0.03
+set terminal pdfcairo color size 7.5 cm, 7.5 cm
+#set terminal pdfcairo color size 15.0 cm, 5.0 cm
+#set multiplot layout 1,4 \
+# margins 0.09,0.97,0.23,0.95 \
+# spacing 0.03,0.03
-unset xlabel
-unset xtics
-set yrange [0:270]
-set ylabel "[kPa]"
-plot \
-"timeseries.txt" u ($0/1008*7):($2/1000) w l lw 2 lc "black" t "effective normal stress", \
-"" u ($0/1008*7):($3/1000) w l lw 2 lt 3 dt "-" t "water pressure"
+#set yrange [4.0:8.0]
-unset xtics
-unset ylabel
-unset ytics
-set yrange [0:4.0]
-set y2range [0:4.0]
-set y2tics 1.0
-set y2label "Shear velocity [km/d]"
-plot "" u ($0/1008*7):($1*60*60*24/1000) w l lw 2 lt 7 t ""
+set key bottom right #samplen 0.9
-set xtics
-unset y2label
-unset y2tics
-set ytics 0.25
-set yrange [1.5:0.0]
-set xlabel "Time [d]"
-set ylabel "Depth of max. strain rate [m]"
-plot "max_strainrate_depth.txt" u ($0/1008*7):(8.0 - $1) w l lw 2 lt 4 t ""
-
-unset multiplot
+set xlabel "Normalized shear velocity [-]"
+set ylabel "Vertical position [m]"
+#set xrange [0.0:1.0]
+#set xtics (0,0.01,0.02)
+#unset ytics
+plot \
+ "two_rivers.txt" u 2:1 w l lw 1 t "simulated Two Rivers till", \
+ "storglaciaren.txt" u 2:1 w l lw 1 t "simulated Storglaciären till", \
+ "upb.txt" u 2:1 w l lw 1 t "simulated Upstream B till"
(DIR) diff --git a/fig3/Makefile b/fig3/Makefile
t@@ -36,10 +36,7 @@ max_strainrate_depth.txt: sim.output00000.txt
"$$f" >> $@; \
done'
-timeseries_combined.txt: timeseries.txt max_strainrate_depth.txt
- paste $^ > $@
-
-../$(FIG).pdf: fig.gp timeseries_combined.txt
+../$(FIG).pdf: fig.gp timeseries.txt max_strainrate_depth.txt
gnuplot $< > $@
clean:
(DIR) diff --git a/fig3/fig.gp b/fig3/fig.gp
t@@ -1,42 +1,36 @@
#!/usr/bin/env gnuplot
reset
-set terminal pdfcairo color size 7.5 cm, 7.5 cm
+set terminal pdfcairo color size 7.5 cm, 8.5 cm
-set multiplot layout 2,1 \
- margins 0.19,0.75,0.23,0.85 \
+set multiplot layout 3,1 \
+ margins 0.19,0.85,0.23,0.95 \
spacing 0.03,0.03
-# set xrange [0.0:0.7]
-# set yrange [0.0:0.01]
-# set key top left #samplen 0.9
-
-set xrange [40:160]
-set x2range [110:-10]
-set x2tics 20 nomirror tc lt 3
+unset xlabel
unset xtics
-
-set x2label "Water pressure [kPa]" tc lt 3
-set ylabel "Shear velocity [km/d]"
-set ytics 1.0
-set arrow from 60,1.0 to 50,2.0 head
-set arrow from 70,3.0 to 83,2.0 head
+set yrange [0:270]
+set ylabel "[kPa]"
plot \
-"timeseries_combined.txt" u ($2/1000):($1*60*60*24/1000) w l lt 7 t ""
+"timeseries.txt" u ($0/1008*7):($2/1000) w l lw 2 lc "black" t "effective normal stress", \
+"" u ($0/1008*7):($3/1000) w l lw 2 lt 3 dt "-" t "water pressure"
-unset arrow
+unset xtics
+unset ylabel
+unset ytics
+set yrange [0:4.0]
+set y2range [0:4.0]
+set y2tics 1.0
+set y2label "Shear velocity [km/d]"
+plot "" u ($0/1008*7):($1*60*60*24/1000) w l lw 2 lt 7 t ""
set xtics
-set xlabel "Effective normal stress [kPa]"
-unset x2tics
-unset ytics
-unset x2label
-unset ylabel
-set y2label "Depth of max. strain rate [m]"
+unset y2label
+unset y2tics
+set ytics 0.25
set yrange [1.5:0.0]
-set y2tics 0.25
-set y2range [1.5:0.0]
-set arrow from 95,0.10 to 85,0.13 head
-set arrow from 100,0.6 to 110,0.8 head
-set arrow from 150,0.8 to 150.3,0.55 head
-plot "" u ($2/1000):(8.0 - $6) w l lt 4 t ""
+set xlabel "Time [d]"
+set ylabel "Depth of max. strain rate [m]"
+plot "max_strainrate_depth.txt" u ($0/1008*7):(8.0 - $1) w l lw 2 lt 4 t ""
+
+unset multiplot
(DIR) diff --git a/fig4/Makefile b/fig4/Makefile
t@@ -3,16 +3,43 @@ FIG = $(shell basename $(CURDIR))
default: ../$(FIG).pdf
-sim.output00000.txt: $(BIN) Makefile ../skindepth.gp
+sim.output00000.txt: $(BIN)
/bin/sh -c '\
./$< --resolution 200 --length 8.0 --normal-stress 150e3 \
--fluid --fluid-permeability 2e-17 --fluid-pressure-top 50e3 \
--fluid-pressure-ampl 50e3 \
--fluid-pressure-freq $$( echo "1.0/(3600*24)" | bc -l ) \
- --file-interval $$( echo "3600" | bc -l ) \
- --time-end $$( echo "3600*24*3 + 3600" | bc -l ) sim'
+ --file-interval $$( echo "60*10" | bc -l ) \
+ --time-end $$( echo "3600*24*7" | bc -l ) sim'
-../$(FIG).pdf: fig.gp sim.output00000.txt
+timeseries.txt: sim.output00000.txt
+ /bin/sh -c '\
+ for f in sim.output*.txt; do \
+ tail -n 1 "$$f" | cut -f2- >> $@; \
+ done'
+
+max_strainrate_depth.txt: sim.output00000.txt
+ /bin/sh -c '\
+ for f in sim.output*.txt; do \
+ awk " \
+ BEGIN{getline; depth=8.0; max=0.0} \
+ NF{ \
+ if (\$$6 >= max && \$$6+0 > 1e-300) \
+ { \
+ max=\$$6; \
+ depth=\$$1; \
+ } \
+ } \
+ END{ \
+ print depth,\"\t\",max \
+ }" \
+ "$$f" >> $@; \
+ done'
+
+timeseries_combined.txt: timeseries.txt max_strainrate_depth.txt
+ paste $^ > $@
+
+../$(FIG).pdf: fig.gp timeseries_combined.txt
gnuplot $< > $@
clean:
(DIR) diff --git a/fig4/fig.gp b/fig4/fig.gp
t@@ -1,62 +1,42 @@
#!/usr/bin/env gnuplot
-
reset
-#set terminal pdfcairo color size 7.5 cm, 7.5 cm
-set terminal pdfcairo color size 15.0 cm, 5.0 cm
-set multiplot layout 1,4 \
- margins 0.09,0.97,0.23,0.95 \
- spacing 0.03,0.03
-
-set yrange [4.0:8.0]
-
-set key bottom right #samplen 0.9
-
-t0 = 47
-tend = 70
+set terminal pdfcairo color size 7.5 cm, 7.5 cm
-filename(t) = sprintf('sim.output%05d.txt', t)
+set multiplot layout 2,1 \
+ margins 0.19,0.75,0.23,0.85 \
+ spacing 0.03,0.03
-set pm3d map
-set palette rgb 10,6,6
+# set xrange [0.0:0.7]
+# set yrange [0.0:0.01]
+# set key top left #samplen 0.9
-load '../skindepth.gp'
-k = 2e-17
-f = 1.0/(3600.0*24.0)
-f(x) = 8.0 - skindepth(f,k)
+set xrange [40:160]
+set x2range [110:-10]
+set x2tics 20 nomirror tc lt 3
+unset xtics
-unset colorbox
-set xlabel "Water pressure [kPa]"
-set ylabel "Vertical position [m]"
-set xrange [0:100]
-plot for [t = t0:tend] filename(t) u ($4/1e3):1:(t-t0) w l lc palette lw 1 t "", \
- f(x) title "" w l
+set x2label "Water pressure [kPa]" tc lt 3
+set ylabel "Shear velocity [km/d]"
+set ytics 1.0
+set arrow from 60,1.0 to 50,2.0 head
+set arrow from 70,3.0 to 83,2.0 head
+plot \
+"timeseries_combined.txt" u ($2/1000):($1*60*60*24/1000) w l lt 7 t ""
-set xlabel "Eff. normal stress [kPa]"
-set ylabel ""
-set xrange [0:200]
-unset ytics
-plot for [t = t0:tend] filename(t) u ($3/1e3):1:(t-t0) w l lc palette lw 1 t "", \
- f(x) title "" w l
+unset arrow
-set xlabel "Shear velocity [m/s]"
-set ylabel ""
-set xrange [0.0:0.06]
-set xtics (0,0.02,0.04,0.06)
+set xtics
+set xlabel "Effective normal stress [kPa]"
+unset x2tics
unset ytics
-plot for [t = t0:tend] filename(t) u 2:1:(t-t0) w l lc palette lw 1 t "", \
- f(x) title "" w l
-
-set xlabel "Shear strain rate [1/s]"
-set ylabel ""
-set xrange [0:3.1e-3]
-#set format x "%.0s*10^{%T}"
-set xtics (0,0.001,0.002,0.003)
-unset ytics
-prev1 = prev2 = 0
-shift(x) = (prev2 = prev1, prev1 = x)
-plot for [t = t0:tend] filename(t) u (shift($2), $0 < 1 ? 1/0 : $2 - prev2):1:(t-t0) w l lc palette lw 1 t "", \
- f(x) title "" w l
-unset format x
-
-unset multiplot
+unset x2label
+unset ylabel
+set y2label "Depth of max. strain rate [m]"
+set yrange [1.5:0.0]
+set y2tics 0.25
+set y2range [1.5:0.0]
+set arrow from 95,0.10 to 85,0.13 head
+set arrow from 100,0.6 to 110,0.8 head
+set arrow from 150,0.8 to 150.3,0.55 head
+plot "" u ($2/1000):(8.0 - $6) w l lt 4 t ""
(DIR) diff --git a/fig5/Makefile b/fig5/Makefile
t@@ -10,8 +10,7 @@ sim.output00000.txt: $(BIN) Makefile ../skindepth.gp
--fluid-pressure-ampl 50e3 \
--fluid-pressure-freq $$( echo "1.0/(3600*24)" | bc -l ) \
--file-interval $$( echo "3600" | bc -l ) \
- --time-end $$( echo "3600*24*3 + 3600" | bc -l ) \
- --normalize sim'
+ --time-end $$( echo "3600*24*3 + 3600" | bc -l ) sim'
../$(FIG).pdf: fig.gp sim.output00000.txt
gnuplot $< > $@
(DIR) diff --git a/fig5/fig.gp b/fig5/fig.gp
t@@ -39,20 +39,19 @@ unset ytics
plot for [t = t0:tend] filename(t) u ($3/1e3):1:(t-t0) w l lc palette lw 1 t "", \
f(x) title "" w l
-set xlabel "Norm. shear vel. [-]"
+set xlabel "Shear velocity [m/s]"
set ylabel ""
-set xrange [0.0:1.0]
-#set xtics (0,0.01,0.02)
+set xrange [0.0:0.06]
+set xtics (0,0.02,0.04,0.06)
unset ytics
plot for [t = t0:tend] filename(t) u 2:1:(t-t0) w l lc palette lw 1 t "", \
f(x) title "" w l
set xlabel "Shear strain rate [1/s]"
set ylabel ""
-#set xrange [0:3.1e-3]
-set xrange [0.0:0.15]
-set xtics (0,0.05,0.10,0.15)
+set xrange [0:3.1e-3]
#set format x "%.0s*10^{%T}"
+set xtics (0,0.001,0.002,0.003)
unset ytics
prev1 = prev2 = 0
shift(x) = (prev2 = prev1, prev1 = x)
(DIR) diff --git a/fig6/Makefile b/fig6/Makefile
t@@ -1,11 +1,23 @@
+BIN = ../1d_fd_simple_shear/1d_fd_simple_shear
FIG = $(shell basename $(CURDIR))
default: ../$(FIG).pdf
-../$(FIG).pdf: fig.gp
+sim.output00000.txt: $(BIN) Makefile ../skindepth.gp
+ /bin/sh -c '\
+ ./$< --resolution 200 --length 8.0 --normal-stress 150e3 \
+ --fluid --fluid-permeability 2e-17 --fluid-pressure-top 50e3 \
+ --fluid-pressure-ampl 50e3 \
+ --fluid-pressure-freq $$( echo "1.0/(3600*24)" | bc -l ) \
+ --file-interval $$( echo "3600" | bc -l ) \
+ --time-end $$( echo "3600*24*3 + 3600" | bc -l ) \
+ --normalize sim'
+
+../$(FIG).pdf: fig.gp sim.output00000.txt
gnuplot $< > $@
clean:
+ $(RM) *.txt
$(RM) ../$(FIG).pdf
.PHONY: default clean
(DIR) diff --git a/fig6/fig.gp b/fig6/fig.gp
t@@ -2,47 +2,62 @@
reset
-set terminal pdfcairo color size 7.5 cm, 6.0 cm
+#set terminal pdfcairo color size 7.5 cm, 7.5 cm
+set terminal pdfcairo color size 15.0 cm, 5.0 cm
+set multiplot layout 1,4 \
+ margins 0.09,0.97,0.23,0.95 \
+ spacing 0.03,0.03
-set lmargin at screen 0.24
-set rmargin at screen 0.75
-set bmargin at screen 0.20
-set tmargin at screen 0.95
+set yrange [4.0:8.0]
set key bottom right #samplen 0.9
-# hydraulic parameters
-phi = 0.25
-beta_f = 3.9e-10
-mu_f = 1.787e-3
-#k = 2e-17
-k_min = 1e-19
-k_max = 1e-13
-
-# forcing parameters
-#f = 1.0/(3600.0*24.0)
-#f_min = 1.0/(3600.0*24*365) # 1/s
-#f_max = 1.0/3600.0 # 1/s
-f_min = 1.0/(24*365) # 1/h
-f_max = 1.0 # 1/h
-
-# f: frequency, k: permeability
-#skindepth(f,k) = (k/(phi*mu_f*beta_f*3.141592654*f))**(0.5) # f: 1/s
-skindepth(f,k) = (k/(phi*mu_f*beta_f*3.141592654*(f/3600)))**(0.5) # f: 1/h
-
-set xlabel "Forcing frequency [1/h]"
-set ylabel "Permeability [m^2]"
-
-set logscale xyzcb
-set xrange [f_min:f_max]
-set yrange [k_min:k_max]
-set samples 100
-set isosamples 60
-
-# Generate contours
-set view map
-unset surface
-set contour
-set cntrparam levels discrete 0.1, 0.3, 1, 3, 10, 30, 100, 300
-set key outside title "Skin depth\n[m]" width 1.0 samplen 2.5
-splot skindepth(x,y) title ""
+t0 = 47
+tend = 70
+
+filename(t) = sprintf('sim.output%05d.txt', t)
+
+set pm3d map
+set palette rgb 10,6,6
+
+load '../skindepth.gp'
+k = 2e-17
+f = 1.0/(3600.0*24.0)
+f(x) = 8.0 - skindepth(f,k)
+
+unset colorbox
+set xlabel "Water pressure [kPa]"
+set ylabel "Vertical position [m]"
+set xrange [0:100]
+plot for [t = t0:tend] filename(t) u ($4/1e3):1:(t-t0) w l lc palette lw 1 t "", \
+ f(x) title "" w l
+
+set xlabel "Eff. normal stress [kPa]"
+set ylabel ""
+set xrange [0:200]
+unset ytics
+plot for [t = t0:tend] filename(t) u ($3/1e3):1:(t-t0) w l lc palette lw 1 t "", \
+ f(x) title "" w l
+
+set xlabel "Norm. shear vel. [-]"
+set ylabel ""
+set xrange [0.0:1.0]
+#set xtics (0,0.01,0.02)
+unset ytics
+plot for [t = t0:tend] filename(t) u 2:1:(t-t0) w l lc palette lw 1 t "", \
+ f(x) title "" w l
+
+set xlabel "Shear strain rate [1/s]"
+set ylabel ""
+#set xrange [0:3.1e-3]
+set xrange [0.0:0.15]
+set xtics (0,0.05,0.10,0.15)
+#set format x "%.0s*10^{%T}"
+unset ytics
+prev1 = prev2 = 0
+shift(x) = (prev2 = prev1, prev1 = x)
+plot for [t = t0:tend] filename(t) u (shift($2), $0 < 1 ? 1/0 : $2 - prev2):1:(t-t0) w l lc palette lw 1 t "", \
+ f(x) title "" w l
+unset format x
+
+unset multiplot
(DIR) diff --git a/fig7/Makefile b/fig7/Makefile
t@@ -1,47 +1,11 @@
-BIN = ../1d_fd_simple_shear/1d_fd_simple_shear
FIG = $(shell basename $(CURDIR))
-RESOLUTION := 1000
-LENGTH := 1.0
-NORMAL_STRESS := 85e3
-
default: ../$(FIG).pdf
-two_rivers.txt: $(BIN) Makefile
- /bin/sh -c '\
- ./$< --resolution $(RESOLUTION) \
- --length $(LENGTH) \
- --normal-stress $(NORMAL_STRESS) \
- --friction-coefficient 0.321 --cohesion 14.0e3 \
- --stress-ratio 0.5 \
- --grain-size 0.1e-3 \
- --normalize' > $@
-
-storglaciaren.txt: $(BIN) Makefile
- /bin/sh -c '\
- ./$< --resolution $(RESOLUTION) \
- --length $(LENGTH) \
- --normal-stress $(NORMAL_STRESS) \
- --friction-coefficient 0.494 --cohesion 5.0e3 \
- --stress-ratio 0.594 \
- --grain-size 0.1e-3 \
- --normalize' > $@
-
-upb.txt: $(BIN) Makefile
- /bin/sh -c '\
- ./$< --resolution $(RESOLUTION) \
- --length $(LENGTH) \
- --normal-stress $(NORMAL_STRESS) \
- --friction-coefficient 0.443 --cohesion 3.0e3 \
- --stress-ratio 0.543 \
- --grain-size 0.1e-3 \
- --normalize' > $@
-
-../$(FIG).pdf: fig.gp two_rivers.txt storglaciaren.txt upb.txt
+../$(FIG).pdf: fig.gp
gnuplot $< > $@
clean:
- $(RM) *.txt
$(RM) ../$(FIG).pdf
.PHONY: default clean
(DIR) diff --git a/fig7/fig.gp b/fig7/fig.gp
t@@ -2,22 +2,47 @@
reset
-set terminal pdfcairo color size 7.5 cm, 7.5 cm
-#set terminal pdfcairo color size 15.0 cm, 5.0 cm
-#set multiplot layout 1,4 \
-# margins 0.09,0.97,0.23,0.95 \
-# spacing 0.03,0.03
+set terminal pdfcairo color size 7.5 cm, 6.0 cm
-#set yrange [4.0:8.0]
+set lmargin at screen 0.24
+set rmargin at screen 0.75
+set bmargin at screen 0.20
+set tmargin at screen 0.95
set key bottom right #samplen 0.9
-set xlabel "Normalized shear velocity [-]"
-set ylabel "Vertical position [m]"
-#set xrange [0.0:1.0]
-#set xtics (0,0.01,0.02)
-#unset ytics
-plot \
- "two_rivers.txt" u 2:1 w l lw 1 t "simulated Two Rivers till", \
- "storglaciaren.txt" u 2:1 w l lw 1 t "simulated Storglaciären till", \
- "upb.txt" u 2:1 w l lw 1 t "simulated Upstream B till"
+# hydraulic parameters
+phi = 0.25
+beta_f = 3.9e-10
+mu_f = 1.787e-3
+#k = 2e-17
+k_min = 1e-19
+k_max = 1e-13
+
+# forcing parameters
+#f = 1.0/(3600.0*24.0)
+#f_min = 1.0/(3600.0*24*365) # 1/s
+#f_max = 1.0/3600.0 # 1/s
+f_min = 1.0/(24*365) # 1/h
+f_max = 1.0 # 1/h
+
+# f: frequency, k: permeability
+#skindepth(f,k) = (k/(phi*mu_f*beta_f*3.141592654*f))**(0.5) # f: 1/s
+skindepth(f,k) = (k/(phi*mu_f*beta_f*3.141592654*(f/3600)))**(0.5) # f: 1/h
+
+set xlabel "Forcing frequency [1/h]"
+set ylabel "Permeability [m^2]"
+
+set logscale xyzcb
+set xrange [f_min:f_max]
+set yrange [k_min:k_max]
+set samples 100
+set isosamples 60
+
+# Generate contours
+set view map
+unset surface
+set contour
+set cntrparam levels discrete 0.1, 0.3, 1, 3, 10, 30, 100, 300
+set key outside title "Skin depth\n[m]" width 1.0 samplen 2.5
+splot skindepth(x,y) title ""