Skip to content

Commit ab4fd58

Browse files
jannitzbonjannitzbon
committed
Initial commit. This version of CryoGrid3 has been used for simulations for an article on ice-wedge degradation.
0 parents  commit ab4fd58

File tree

89 files changed

+3776
-0
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

89 files changed

+3776
-0
lines changed

‎CryoGrid3_xice_mpi.m‎

Lines changed: 369 additions & 0 deletions
Large diffs are not rendered by default.

‎README.md‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
# CryoGrid3
2+
CryoGrid 3 is a simple land-surface scheme dedicated to modeling of ground temperatures in permafrost environments.
3+
4+
The model version of is the refernce for the research article "Modelling the degradation of ice-wedges in polygonal tundra under different hydrological conditions" which is intended to be published in the scientific journal "The Cryosphere".
5+
6+
The parameters are set to the default values used in the article. Parameters different from the default values can be specified in the main scirpt "CryoGrid3_xice_mpi" (general parameters) and in the function "get_parallel_variables" (tile-specific parameters) which is part of the module "cryoGridLateral".
7+
8+
To start the program, just run the script "CryoGrid3_xice_mpi". The default output directory is "./runs/".

‎add_modules.m‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,13 @@
1+
clear all
2+
close all
3+
4+
% import CryoGrid modules
5+
addpath('modules/cryoGridTechnical/')
6+
addpath('modules/cryoGridInitialize/')
7+
addpath('modules/cryoGridSEB/')
8+
addpath('modules/cryoGridSoil/')
9+
addpath('modules/cryoGridSnow/')
10+
addpath('modules/cryoGridInfiltrationUnfrozenSoil')
11+
addpath('modules/cryoGridExcessIce/')
12+
addpath('modules/cryoGridExcessIceInfiltration')
13+
addpath('modules/cryoGridLateral')

‎forcing/samoylov_ERA_obs_fitted_1979_2014_spinup_extended2044.mat‎

17.3 MB
Binary file not shown.

‎modules/cryoGridExcessIce/excessGroundIce.m‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
function [GRID, PARA] = excessGroundIce(T, GRID, PARA)
2+
3+
if ~isempty(PARA.soil.mobileWaterDomain) && (sum(double(T(GRID.soil.cT_domain)>0 & GRID.soil.excessGroundIce==1))~=0) && isempty(GRID.snow.cT_domain_ub)
4+
disp('excess ice thawing');
5+
GRID.soil.excessGroundIce = GRID.soil.excessGroundIce==1 & T(GRID.soil.cT_domain)<=0; %remove the thawed cell from the list
6+
[GRID meltwaterGroundIce PARA] = excessGroundIceThaw4(T, GRID, PARA); %meltwaterGroundIce could be read out, but is not yet implemented
7+
GRID = updateGRID_excessice(GRID);
8+
end

‎modules/cryoGridExcessIce/excessGroundIceThaw4.m‎

Lines changed: 140 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,140 @@
1+
function [GRID, meltwaterGroundIce, PARA]=excessGroundIceThaw4(T, GRID, PARA)
2+
3+
%disp('rearranging grid cells due to ground ice thaw')
4+
5+
waterLevel=PARA.soil.waterTable; %remove supersaturation only when there is no snow on top of soil!!!!!!!
6+
7+
8+
mineral=GRID.general.K_delta(GRID.soil.cT_domain).*GRID.soil.cT_mineral; %calculates amounts in [m]
9+
organic=GRID.general.K_delta(GRID.soil.cT_domain).*GRID.soil.cT_organic;
10+
water=GRID.general.K_delta(GRID.soil.cT_domain).*GRID.soil.cT_water;
11+
natPor=GRID.general.K_delta(GRID.soil.cT_domain).*GRID.soil.cT_natPor;
12+
13+
cT_grid=GRID.general.cT_grid(GRID.soil.cT_domain);
14+
K_delta=GRID.general.K_delta(GRID.soil.cT_domain);
15+
16+
17+
mobileWater = double(T(GRID.soil.cT_domain)>0) .* (water-natPor) .* double(water>natPor);
18+
[startCell ~]= LayerIndex(mobileWater~=0); %this is faster
19+
20+
%move solids down
21+
for i=startCell:-1:1
22+
F_solid_down=K_delta(i)-mineral(i)-organic(i)-natPor(i);
23+
j=i-1;
24+
while j>0 && F_solid_down>0
25+
mineralDown = min(mineral(j), mineral(j)./(mineral(j)+organic(j)).*F_solid_down);
26+
organicDown = min(organic(j), organic(j)./(mineral(j)+organic(j)).*F_solid_down);
27+
mineral(i)=mineral(i)+mineralDown;
28+
organic(i)=organic(i)+organicDown;
29+
mineral(j)=mineral(j)-mineralDown;
30+
organic(j)=organic(j)-organicDown;
31+
F_solid_down=F_solid_down-mineralDown-organicDown;
32+
j=j-1;
33+
end
34+
end
35+
36+
%adjust the natural porosity
37+
natPor(1:startCell)=K_delta(1:startCell)-mineral(1:startCell)-organic(1:startCell);
38+
39+
%move water up
40+
mobileWater=0;
41+
for i=startCell:-1:1
42+
totalWater=water(i)+mobileWater;
43+
mobileWater=totalWater-natPor(i);
44+
mobileWater=max(0,mobileWater);
45+
water(i)=totalWater-mobileWater;
46+
end
47+
48+
%clean up grid cells with non-zero+non-unity water content in domains without soil matrix
49+
mobileWater=0;
50+
for i=1:startCell
51+
if mineral(i)+organic(i)==0
52+
mobileWater=mobileWater+water(i);
53+
water(i)=0;
54+
end
55+
end
56+
for i=startCell:-1:1
57+
if mineral(i)+organic(i)==0
58+
water(i)=min(K_delta(i), mobileWater);
59+
mobileWater=mobileWater-water(i);
60+
water(i)=round(water(i)./K_delta(i)).*K_delta(i); %this violates the water balance, but ensures that no grid cells with partly water and partly air can exist;
61+
end
62+
end
63+
64+
65+
GRID.soil.cT_mineral=mineral./K_delta;
66+
GRID.soil.cT_organic=organic./K_delta;
67+
GRID.soil.cT_water=water./K_delta;
68+
GRID.soil.cT_natPor=natPor./K_delta;
69+
test=mineral+water+organic;
70+
GRID.soil.cT_soilType(water(:,1)==1,1)=1; %sets sand freeze curve for all water grid cells (comment: find was replaced!!!)
71+
72+
GRID.soil.K_mineral(1)=GRID.soil.cT_mineral(1);
73+
GRID.soil.K_mineral(2:startCell+1)=(GRID.soil.cT_mineral(2:startCell+1)+GRID.soil.cT_mineral(1:startCell))/2 ;
74+
GRID.soil.K_organic(1)=GRID.soil.cT_organic(1);
75+
GRID.soil.K_organic(2:startCell+1)=(GRID.soil.cT_organic(2:startCell+1)+GRID.soil.cT_organic(1:startCell))/2 ;
76+
GRID.soil.K_water(1)=GRID.soil.cT_water(1);
77+
GRID.soil.K_water(2:startCell+1)=(GRID.soil.cT_water(2:startCell+1)+GRID.soil.cT_water(1:startCell))/2 ;
78+
% GRID.soil.K_natPor(1)=GRID.soil.cT_natPor(1);
79+
% GRID.soil.K_natPor(2:startCell+1)=(GRID.soil.cT_natPor(2:startCell+1)+GRID.soil.cT_natPor(1:startCell))/2 ;
80+
GRID.soil.K_soilType(1)=GRID.soil.cT_soilType(1);
81+
GRID.soil.K_soilType(2:startCell+1)=round((GRID.soil.cT_soilType(2:startCell+1)+GRID.soil.cT_soilType(1:startCell))/2) ;
82+
83+
meltwaterGroundIce=0;
84+
85+
while (GRID.soil.cT_mineral(1)+GRID.soil.cT_organic(1)+GRID.soil.cT_water(1)==0) || (GRID.soil.cT_water(1)==1 && GRID.general.K_grid(GRID.soil.K_domain_ub)<waterLevel) %remove grid cells until the water level is reached (comment: why GRID.general.K_grid(GRID.soil.K_domain_ub)<waterLevel)
86+
87+
GRID.air.cT_domain(GRID.soil.cT_domain_ub)=1;
88+
GRID.air.K_domain(GRID.soil.K_domain_ub)=1;
89+
GRID.air.cT_domain_lb=GRID.air.cT_domain_lb+1;
90+
GRID.air.K_domain_lb=GRID.air.K_domain_lb+1;
91+
92+
GRID.soil.cT_domain(GRID.soil.cT_domain_ub)=0;
93+
GRID.soil.K_domain(GRID.soil.K_domain_ub)=0;
94+
GRID.soil.soilGrid(1)=[];
95+
%GRID.soil.convectiveDomain(1)=[];
96+
GRID.soil.cT_water(1)=[];
97+
GRID.soil.cT_organic(1)=[];
98+
GRID.soil.cT_natPor(1)=[];
99+
GRID.soil.cT_mineral(1)=[];
100+
GRID.soil.cT_soilType(1)=[];
101+
GRID.soil.K_water(1)=[];
102+
GRID.soil.K_organic(1)=[];
103+
GRID.soil.K_mineral(1)=[];
104+
GRID.soil.K_soilType(1)=[];
105+
GRID.soil.excessGroundIce(1)=[];
106+
meltwaterGroundIce=meltwaterGroundIce+GRID.general.K_delta(GRID.soil.cT_domain_ub);
107+
GRID.soil.cT_domain_ub=GRID.soil.cT_domain_ub+1;
108+
GRID.soil.K_domain_ub=GRID.soil.K_domain_ub+1;
109+
end
110+
111+
GRID = initializeSoilThermalProperties( GRID, PARA );
112+
113+
%reduce water content above the perched water table
114+
%it might be also good to use an external function for this since this is
115+
%only an option in the model
116+
%i=0;
117+
% while i./startCell < 1-PARA.soil.perchedWaterTable
118+
% GRID.soil.cT_water(i+1,1) = PARA.soil.saturationAbovePerchedWaterTable.*(1-GRID.soil.cT_mineral(i+1,1)-GRID.soil.cT_organic(i+1,1));
119+
%
120+
% GRID.soil.K_water(i+1,1) = PARA.soil.saturationAbovePerchedWaterTable.*(1-GRID.soil.K_mineral(i+1,1)-GRID.soil.K_organic(i+1,1));
121+
% i=i+1;
122+
% end
123+
124+
125+
% [GRID.soil.cT_frozen,...
126+
% GRID.soil.cT_thawed,...
127+
% GRID.soil.K_frozen,...
128+
% GRID.soil.K_thawed,...
129+
% GRID.soil.conductivity,...
130+
% GRID.soil.capacity] = initialize(GRID.soil.cT_water,...
131+
% GRID.soil.cT_mineral,...
132+
% GRID.soil.cT_organic,...
133+
% GRID.soil.cT_soilType,...
134+
% GRID.soil.K_water,...
135+
% GRID.soil.K_mineral,...
136+
% GRID.soil.K_organic,...
137+
% GRID.soil.K_soilType,...
138+
% PARA.technical.arraySizeT,...
139+
% GRID.general.cT_grid(GRID.soil.cT_domain),...
140+
% PARA.soil.kh_bedrock);

‎modules/cryoGridExcessIce/initializeExcessIce2.m‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
function [GRID,PARA] = initializeExcessIce2(GRID,PARA)
2+
3+
GRID.soil.excessGroundIce = GRID.soil.cT_water>GRID.soil.cT_natPor;

‎modules/cryoGridExcessIce/mixingWaterBody.m‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
function [T] = mixingWaterBody(T, GRID)
2+
3+
% mixing of temperatures in unfrozen water domain
4+
if GRID.soil.cT_organic(1)+GRID.soil.cT_mineral(1)<=1e-6 && T(GRID.soil.cT_domain_ub)>0
5+
unfrozenWaterBody = GRID.soil.cT_organic+GRID.soil.cT_mineral<=1e-6 & T(GRID.soil.cT_domain)>0;
6+
waterCells = sum( unfrozenWaterBody );
7+
if waterCells > 1
8+
Tav = sum( T(GRID.soil.cT_domain_ub:GRID.soil.cT_domain_ub+waterCells-1) .* GRID.general.K_delta(GRID.soil.cT_domain_ub:GRID.soil.cT_domain_ub+waterCells-1) ) ...
9+
./ sum( GRID.general.K_delta(GRID.soil.cT_domain_ub:GRID.soil.cT_domain_ub+waterCells-1) ) ;
10+
T(GRID.soil.cT_domain_ub:GRID.soil.cT_domain_ub+waterCells-1) = Tav;
11+
end
12+
end
13+
14+
end

‎modules/cryoGridExcessIce/updateGRID_excessice.m‎

Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,52 @@
1+
function GRID = updateGRID_excessice(GRID)
2+
%--- update soil and water grid after subsidence ----------------------
3+
4+
% change soil with 100% water to water cell
5+
soilGRIDsize = sum(GRID.soil.cT_domain);
6+
7+
% set all cells above without soil to air (water runs off and no snow cover
8+
% present)
9+
GRID.air.cT_domain(GRID.soil.cT_domain) = (GRID.soil.cT_organic==0 & GRID.soil.cT_mineral==0);
10+
GRID.air.K_domain(GRID.soil.K_domain) = (GRID.soil.K_organic==0 & GRID.soil.K_mineral==0);
11+
12+
GRID.soil.cT_domain(GRID.soil.cT_domain) = (GRID.soil.cT_organic>0 | GRID.soil.cT_mineral>0);
13+
GRID.soil.K_domain(GRID.soil.K_domain) = (GRID.soil.K_organic>0 | GRID.soil.K_mineral>0);
14+
15+
if soilGRIDsize ~= sum(GRID.soil.cT_domain)
16+
17+
disp('subsidence - updating grid information');
18+
19+
%water_bucket = GRID.soil.cT_water(1);
20+
cT_no_water = (GRID.soil.cT_organic>0 | GRID.soil.cT_mineral>0);
21+
K_no_water = (GRID.soil.K_organic>0 | GRID.soil.K_mineral>0);
22+
23+
[GRID.soil.cT_domain_lb GRID.soil.cT_domain_ub] = LayerIndex(GRID.soil.cT_domain);
24+
[GRID.soil.K_domain_lb GRID.soil.K_domain_ub] = LayerIndex(GRID.soil.K_domain);
25+
26+
[GRID.air.cT_domain_lb GRID.air.cT_domain_ub] = LayerIndex(GRID.air.cT_domain);
27+
[GRID.air.K_domain_lb GRID.air.K_domain_ub] = LayerIndex(GRID.air.K_domain);
28+
29+
%-- update all other soil grid infos if size has changed
30+
% adjust cT grid fields
31+
GRID.soil.cT_water = GRID.soil.cT_water(cT_no_water);
32+
GRID.soil.cT_mineral = GRID.soil.cT_mineral(cT_no_water);
33+
GRID.soil.cT_organic = GRID.soil.cT_organic(cT_no_water);
34+
GRID.soil.cT_soilType = GRID.soil.cT_soilType(cT_no_water);
35+
GRID.soil.cT_natPor = GRID.soil.cT_natPor(cT_no_water);
36+
GRID.soil.excessGroundIce = GRID.soil.excessGroundIce(cT_no_water);
37+
GRID.soil.conductivity = GRID.soil.conductivity(cT_no_water, :);
38+
GRID.soil.capacity = GRID.soil.capacity(cT_no_water, :);
39+
GRID.soil.liquidWaterContent = GRID.soil.liquidWaterContent(cT_no_water, :);
40+
GRID.soil.cT_frozen = GRID.soil.cT_frozen(cT_no_water);
41+
GRID.soil.cT_thawed = GRID.soil.cT_thawed(cT_no_water);
42+
GRID.soil.K_frozen = GRID.soil.cT_frozen; % this is ok since K_frozen and cT_frozen
43+
GRID.soil.K_thawed = GRID.soil.cT_thawed; % are the same from the start (see initialize.m)
44+
45+
% adjust K grid fields
46+
GRID.soil.soilGrid = GRID.soil.soilGrid(K_no_water);
47+
GRID.soil.K_water = GRID.soil.K_water(K_no_water);
48+
GRID.soil.K_mineral = GRID.soil.K_mineral(K_no_water);
49+
GRID.soil.K_organic = GRID.soil.K_organic(K_no_water);
50+
GRID.soil.K_soilType = GRID.soil.K_soilType(K_no_water);
51+
52+
end

‎modules/cryoGridExcessIceInfiltration/excessGroundIceInfiltration.m‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
function [GRID, PARA, wc, meltwaterGroundIce] = excessGroundIceInfiltration(T, wc, GRID, PARA)
2+
meltwaterGroundIce=0;
3+
if ~isempty(PARA.soil.mobileWaterDomain) && (sum(double(T(GRID.soil.cT_domain)>0 & GRID.soil.excessGroundIce==1))~=0) && isempty(GRID.snow.cT_domain_ub)
4+
disp('excessGroundIceInfiltration - excess ice thawing');
5+
GRID.soil.excessGroundIce = GRID.soil.excessGroundIce==1 & T(GRID.soil.cT_domain)<=0; %remove the thawed cell from the list
6+
[GRID, meltwaterGroundIce, wc] = excessGroundIceThaw4Infiltration(T, wc, GRID, PARA); %meltwaterGroundIce could be read out, but is not yet implemented
7+
8+
end

0 commit comments

Comments
 (0)