Skip to content

Commit

Permalink
Initial commit version from 15th Feb 2021
Browse files Browse the repository at this point in the history
  • Loading branch information
@NiklasHohmann
NiklasHohmann authored Mar 16, 2023
0 parents commit 9505c6e
Show file tree
Hide file tree
Showing 54 changed files with 6,772 additions and 0 deletions.
13 changes: 13 additions & 0 deletions README.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
double click on the file carboCATGUI.m to open Matlab and load this file.
Run the function either by enterting the function name carboCATGUI at the command prompt,
or by clicking run in the editor tab

The GUI should be reasonably intuitive; click on the Initialise button to set up the model
ready to run, then click on Run CA Model, then on PLot run and you should see the
results from the paramsProcess.txt and the paramsInputValues.txt input files I have set up
in the params folder. You can edit both of these files to run different models.
You can also use the code in the utilities folder to make new initial topography,
intial facies maps, sea-level curves and other inputs you might need to create
useful models

Any questions, e-mail peter.burgess@liverpool.ac.uk
207 changes: 207 additions & 0 deletions calculateErrorBetweenMarginTrajectory.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,207 @@
clear all
clc

%%load saved data
strata=load('C:\EquinorProject\Meeting\28Apr\WaveNorth\strata.mat');
glob.strata=strata.strata;
numberOfLayers=load('C:\EquinorProject\Meeting\28Apr\WaveNorth\numberOfLayers.mat');
glob.numberOfLayers=numberOfLayers.numberOfLayers;
facies=load('C:\EquinorProject\Meeting\28Apr\WaveNorth\facies.mat');
glob.facies=facies.facies;
glob.ySize=size(glob.strata,1);
glob.xSize=size(glob.strata,2);
iteration=size(glob.strata,3);
glob.dx=250; %cell size
deltaT=100; %calculate platform margin at every 100 iterations
T=linspace(deltaT,iteration,iteration/deltaT);

% %%for diagonal cross section where y=x;
% long=min(glob.xSize,glob.ySize);
% gradDiag=zeros(long,iteration);
% PfDiag=zeros(long,iteration);
% Zdiag=zeros(long,iteration);
% for y=1:long-1
% xPos=y;
% for t=2:iteration
%
% g=glob.strata(y+1,xPos+1,t)-glob.strata(y,xPos,t);
% gradDiag(y,t)=g/25;
% for m=1:glob.numberOfLayers(y,xPos,t)
% if glob.facies{y+1,xPos,t}(1)==10 && glob.facies{y,xPos,t}(m)==1
% PfDiag(y,t)=1;
% Zdiag(y,t)=glob.strata(y,xPos,t);
% elseif glob.facies{y,xPos,t}(m)==10 && glob.facies{y,xPos,t}(m)==11 && glob.facies{y+1,xPos,t}(1)==1
% PfDiag(y,t)=1;
% end
% end
%
%
% end
%
% end
% gradDiag(gradDiag>1)=1;
% gradDiag=abs(gradDiag);
% for i=1:length(T)
% t=T(i);
% PDiag(:,i)=gradDiag(:,t).*PfDiag(:,t);
% end
% [Xdiag,Zdiag] = find(PDiag);
% for i=1:size(Xdiag,1)
% Z2diag(i)=glob.strata(Xdiag(i),Xdiag(i),Zdiag(i)*deltaT);
% end
% figure
% subplot(2,1,1)
% scatter(Xdiag,Zdiag,200,'r');
% ylabel('EMT','FontSize',18,'FontWeight','bold');
% xlabel('km','FontSize',18,'FontWeight','bold');
% title(['Diagonal, t= ', num2str(0.001*iteration),'My'],'FontSize',20);
% xticks([1, 0.25*(long+2), 0.5*(long+2), 0.75*(long+2),long+1]);
% xticklabels (round([0 ,long*1.414*glob.dx/4000 ,long*1.414*glob.dx/2000,long*1.414*3*glob.dx/4000,long*1.414*glob.dx/1000]));
% % yticks([0 ,5, 10]);
% % yticklabels ([0 0.5*iteration/1000,iteration/1000]);
% axis([1 long+1 0 iteration/deltaT])
% ax=gca;
% ax.LineWidth = 0.6;
% ax.FontSize = 24;
% ax.FontWeight = 'bold';
%
% subplot(2,1,2)
% scatter(Xdiag,Z2diag,200,'r');
% ylabel('Elevation (m)','FontSize',18,'FontWeight','bold');
% xlabel('km','FontSize',18,'FontWeight','bold');
% xticks([1, 0.25*(long+2), 0.5*(long+2), 0.75*(long+2),long+1]);
% xticklabels (round([0 ,long*1.414*glob.dx/4000 ,long*1.414*glob.dx/2000,long*1.414*3*glob.dx/4000,long*1.414*glob.dx/1000]));
% % yticks([0 ,5, 10]);
% % yticklabels ([0 0.5*iteration/1000,iteration/1000]);
% axis([1 long+1 0 1500])
% ax=gca;
% ax.LineWidth = 0.6;
% ax.FontSize = 24;
% ax.FontWeight = 'bold';


%%for othogonal cross section along y axis
xPos=glob.xSize/2;
gradY=zeros(glob.ySize,iteration); %to store elevation difference between adjacent cells
Pfy=zeros(glob.ySize,iteration); %to store probability of facies transition
Zy=zeros(glob.ySize,iteration); %to store elevation values of platform margins
for y=1:glob.ySize-1

for t=2:iteration

g=glob.strata(y+1,xPos,t)-glob.strata(y,xPos,t);
gradY(y,t)=g/50;
for m=1:glob.numberOfLayers(y,xPos,t)
if glob.facies{y+1,xPos,t}(1)==10 && glob.facies{y,xPos,t}(m)==1
Pfy(y,t)=1;
Zy(y,t)=glob.strata(y,xPos,t);
elseif glob.facies{y,xPos,t}(m)==10 && glob.facies{y,xPos,t}(m)==11 && glob.facies{y+1,xPos,t}(1)==1
Pfy(y,t)=1;
end
end


end
end
%turn the calculated gradient into a value between 0 and 1
gradY(gradY>1)=1;
gradY=abs(gradY);
for i=1:length(T)
t=T(i);
Py(:,i)=gradY(:,t).*Pfy(:,t);
end
[X1,Z] = find(Py);
for i=1:size(X1,1)
Z2(i)=glob.strata(X1(i),xPos,Z(i)*deltaT);
end

%%plot figure in time
figure
subplot('Position',[0.14, 0.55, 0.72 0.35])
scatter(X1,Z,200,'r');
ylabel('EMT','FontSize',18,'FontWeight','bold');
xlabel('y(km)','FontSize',18,'FontWeight','bold');
title(['t= ', num2str(0.001*iteration),'My'],'FontSize',20);
xticks([1 10 20 30 40 50 60 70 80 90 100]);
xticklabels (0:(glob.dx/100):20);
% yticks([0 ,5, 10]);
% yticklabels ([0 0.5*iteration/1000,iteration/1000]);
axis([1 glob.ySize 0 iteration/deltaT])
ax=gca;
ax.LineWidth = 0.6;
ax.FontSize = 24;
ax.FontWeight = 'bold';

%%for othogonal cross section along x axis
yPos=glob.ySize/2;
gradX=zeros(glob.xSize,iteration);
Pfx=zeros(glob.xSize,iteration);
for x=1:glob.xSize-1

for t=2:iteration
g=glob.strata(yPos,x+1,t)-glob.strata(yPos,x,t);
gradX(x,t)=g/50;
for m=1:glob.numberOfLayers(yPos,x,t)
if glob.facies{yPos,x+1,t}(1)==10 && glob.facies{yPos,x,t}(m)==1
Pfx(x,t)=1;
elseif glob.facies{yPos,x,t}(m)==10 && glob.facies{yPos,x,t}(m)==11 && glob.facies{yPos,x+1,t}(1)==1
Pfx(x,t)=1;
end
end


end
end
gradX(gradX>1)=1;
gradX=abs(gradX);
[row, col] = find(Pfx(:,1:iteration)>0);
for i=1:length(T)
t=T(i);
Px(:,i)=gradX(:,t).*Pfx(:,t);
end

[X2,Z] = find(Px);
for i=1:size(X2,1)
Z3(i)=glob.strata(yPos,X2(i),Z(i)*deltaT);
end
subplot('Position',[0.14, 0.1, 0.68 0.35])
scatter(X2,Z,200,'r');
ylabel('EMT','FontSize',18,'FontWeight','bold');
xlabel('x(km)','FontSize',18,'FontWeight','bold');
xticks([1 10 20 30 40 50 60 70 80 90 100]);
xticklabels (0:(glob.dx/100):20);
% yticks([0 ,5, 10]);
% yticklabels ([0 0.5*iteration/1000,iteration/1000]);
axis([1 glob.xSize 0 iteration/deltaT])
ax=gca;
ax.LineWidth = 0.6;
ax.FontSize = 24;
ax.FontWeight = 'bold';


%%plot cross section in depth
figure
subplot('Position',[0.14, 0.55, 0.72 0.35])
scatter(X1,Z2,200,'r');
ylabel('Elevation (m)','FontSize',18,'FontWeight','bold');
xlabel('y(km)','FontSize',18,'FontWeight','bold');
title(['t= ', num2str(0.001*iteration),'My'],'FontSize',20);
xticks([1 10 20 30 40 50 60 70 80 90 100]);
xticklabels (0:(glob.dx/100):20);
axis([1 glob.ySize 0 1500])
ax=gca;
ax.LineWidth = 0.6;
ax.FontSize = 24;
ax.FontWeight = 'bold';

subplot('Position',[0.14, 0.1, 0.68 0.35])
scatter(X2,Z3,200,'r');
ylabel('Elevation (m)','FontSize',18,'FontWeight','bold');
xlabel('x(km)','FontSize',18,'FontWeight','bold');
xticks([1 10 20 30 40 50 60 70 80 90 100]);
xticklabels (0:(glob.dx/100):20);
axis([1 glob.xSize 0 1500])
ax=gca;
ax.LineWidth = 0.6;
ax.FontSize = 24;
ax.FontWeight = 'bold';
174 changes: 174 additions & 0 deletions calculateFaciesCAProdDependent.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,174 @@
function [glob] = calculateFaciesCAProdDependent(glob, iteration)
% Calculate the facies cellular automata according to neighbour rules in
% glob.CARules and modify production rate modifier for each cell according to number of neighbours.
% The order in which facies are checked depends on the production value at the depth of the cell

j = iteration - 1;
k = iteration;

% Neighbours contains a count of neighbours for each facies at each grid point and is populated by the countAllNeighbours function
neighbours = zeros(glob.ySize, glob.xSize, glob.maxProdFacies);

[neighbours] = countAllNeighbours(glob, j, neighbours);
glob.faciesProdAdjust = zeros(glob.ySize, glob.xSize);

tempFaciesProd = zeros(glob.ySize, glob.xSize, glob.maxProdFacies); % need a CA array for each prod facies

for y = 1 : glob.ySize
for x= 1 : glob.xSize

oneFacies = glob.facies{y,x,j}(1);

% Only do anything here if the latest stratal surface is below sea-level, i.e.
% water depth > 0.001 and if the concentration is enough to maintain production
if glob.wd(y,x,k) > 0.001

% For a subaerial hiatus, now reflooded because from above wd > 0
if oneFacies == 9 % 9 is the code for subaerial exposure
checkProcess=strcmp(glob.refloodingRoutine,'pre-exposed');
if checkProcess==1
% Get the facies from the pre-exposed
glob.facies{y,x,k}(1) = findPreHiatusFacies(glob, x,y,iteration); % reoccupy with facies from below hiatus
glob.numberOfLayers(y,x,k)=1;
else
% Set facies =0;
glob.facies{y,x,k}(1) = 0;
glob.numberOfLayers(y,x,k)=0;
end
if glob.facies{y,x,k}(1)> glob.maxProdFacies || glob.facies{y,x,k}(1)==0
glob.facies{y,x,k}(1) = 0;
glob.numberOfLayers(y,x,k)=0;
end

% For cells already containing producing facies
elseif (oneFacies <= glob.maxProdFacies && oneFacies>0) && glob.numberOfLayers(y,x,j)==1

% Check if neighbours is less than min for survival CARules(i,2) or
% greater than max for survival CARules(i,3), or if water depth is greater than production cutoff and if so kill
% that facies
thick = glob.carbonateVolMap(y,x) / (glob.ySize*glob.xSize);
checkProcess=strcmp(glob.concentrationRoutine,'on');

if (neighbours(y,x,oneFacies) < glob.CARules(oneFacies,2)) ||...
(neighbours(y,x,oneFacies) > glob.CARules(oneFacies,3)) ||...
(glob.wd(y,x,iteration) >= glob.prodRateWDCutOff(oneFacies)) ||...
( checkProcess == 1 && thick < 0.01)
glob.facies{y,x,k}(1) = 0;
else % The right number of neighbours exists
glob.facies{y,x,k}(1) = oneFacies;
glob.numberOfLayers(y,x,k)=1;
end

else % Otherwise cell must be empty or contain transported product so see if it can be colonised with a producing facies
% Check again the carbonate concentration
thick = glob.carbonateVolMap(y,x) / (glob.ySize*glob.xSize);
checkProcess=strcmp(glob.concentrationRoutine,'on');
if checkProcess == 0 || thick > 0.01
for f = 1:glob.maxProdFacies
% Check if number of neighbours is within range to trigger new
% facies cell, and only allow new cell if the water depth is less the production cut off depth
if (neighbours(y,x,f) >= glob.CARules(f,4)) && (neighbours(y,x,f) <= glob.CARules(f,5)) && (glob.wd(y,x,iteration) < glob.prodRateWDCutOff(f))
tempFaciesProd(y,x,f) = f; % new facies cell triggered at y,x
end
end

prodRate = zeros(1,glob.maxProdFacies);
for f = 1:glob.maxProdFacies
if tempFaciesProd(y,x,f) > 0
prodRate(f) = glob.prodRate(f) * calculateProductionWDAdjustment(glob, x,y, f, iteration);
end
end

% Find the highest production rate and assign a fraction of 1 at the correct position in the glob.fraction array
[sortedArray,index]=sort(prodRate);
if sum(sortedArray(:))==0
glob.facies{y,x,k}(1)=0;
else
glob.facies{y,x,k}(1)=index(glob.maxProdFacies); % the last element in the array has the largest production
glob.numberOfLayers(y,x,k)=1;
end
end
end

% Finally, calculate the production adjustment factor for the current facies distribution
oneFacies = glob.facies{y,x,k}(1);
glob.faciesProdAdjust(y,x)=calculateFaciesProductionAdjustementFactor(glob,y,x,oneFacies,neighbours);
else
glob.facies{y,x,k}(1) = 9; % Set to above sea-level facies because must be at or above sea-level
glob.numberOfLayers(y,x,k)=1;

end
end
end
end

function [neighbours] = countAllNeighbours(glob, j, neighbours)
% Function to count the number of cells within radius containing facies 1
% to maxFacies across the whole grid and store results in neihgbours matrix

ySize=int16(glob.ySize);
xSize=int16(glob.xSize);
for yco = 1 : ySize
for xco = 1 : xSize
oneFacies = glob.facies{yco,xco,j}(1);
if oneFacies>0 && oneFacies<=glob.maxProdFacies
radius = glob.CARules(oneFacies,1);
else
radius = glob.CARules(1,1);
end
for l = -radius : radius
for m = -radius : radius

y = yco + l;
x = xco + m;

checkProcess=strcmp(glob.wrapRoutine,'unwrap');
if checkProcess==1
%if near the edge, complete in a mirror-like image the
%neighbours array
if y<1; y=1+(1-y); end
if x<1; x=1+(1-x); end
if y>ySize; y=ySize+(ySize-y); end
if x>xSize; x=xSize+(xSize-x); end
else
%or wrap around the corners
if y<1; y=ySize+1+l; end
if x<1; x=xSize+1+m; end
if y>ySize; y=l; end
if x>xSize; x=m; end
end
%now count the neighbours using the x-y indeces

%Don't include cell that has over a certain wd
%difference, and penalise the CA with a factor
wdDiff=abs(glob.wd(yco,xco,j)-glob.wd(y,x,j));

if wdDiff>glob.BathiLimit
wdDiffFactor=0;
else
wdDiffFactor=(-1/glob.BathiLimit)*wdDiff+1;
end
if wdDiffFactor>0.9; wdDiffFactor=1; end
faciesType=glob.facies{y,x,j}(1);
% Count producing facies as neighbours but do not include the center cell -
% neighbours count should not include itself
if faciesType > 0 && faciesType <= glob.maxProdFacies && not (l == 0 && m == 0)
neighbours(yco,xco,faciesType) = neighbours(yco,xco,faciesType) + (1*wdDiffFactor);
end

end
end
end
end
end

function [preHiatusFacies] = findPreHiatusFacies(glob, x,y,iteration)

k = iteration - 1;

while k > 0 && glob.facies{y,x,k}(1) == 9
k = k - 1;
end

preHiatusFacies = glob.facies{y,x,k}(1);
end
Loading

0 comments on commit 9505c6e

Please sign in to comment.