% This Version Modified and Commented on January 24, 2015

% this version modified and commented on January 24, 2015

% set constants

WHC=20; % mm water in holding capacity, in assumed 50 cm depth

PotET=10; % potential ET in time step

MaxCaconc1=1; % Ca concentration that fully inhibits weathering of mineral 1, reflecting solubility control of mineral

% weathering; the three minerals set to 1, .5. .2 for basalt and 1, .1, and

% 1, .1, and .01 for non-basalt

MaxCaconc2=.5; % Ca concentration that fully inhibits weathering of mineral 2

MaxCaconc3=.2; % Ca concentration that fully inhibits weathering of mineral 3

weathercoeff1=1;% kinetic coefficients; also 1, .5. and .2 for basalt and 1, .1, and .01 for non-basalt.

weathercoeff2=0.5;

weathercoeff3=.2;

variablerain=1; %if this is set to 1, rainfall incorporates random variation, 0 it is determinate

begin=2;

stop=100000; %number of time steps for simulation

startaverage=99000; %average results for last interval; important when there is a random component to rainfall

changerain=10000001; %if you want to change rainfall within a run, to explore reversibility, this is the year it changes;

% if you don't want to change set it outside the run

changeamnt=.5; %the factor by which rainfall changes; can increase or decrease

% initialize arrays for outer loop of rainfall

waterm(1)=1; % these are all terms that are averaged for the last interval of a run

rain(1)=1;

meanCaconc(1)=1;

meanproportiontotalCarem(1)=1;

meanrain(1)=1;

meanweather(1)=1;

year(1)=1; % The initial year of simulation

%Here set the first element of the output matrix for the properties

%you may want to graph, equal to the value in year 1 above. I use this

%when outputting results to excel, when I'm going to use the output within another package.

graph(1,1)=waterm(1);

graph(2,1)=meanrain(1);

graph(3,1)=meanCaconc(1);

graph(4,1)=meanproportiontotalCarem(1);

graph(5,1)=meanweather(1);

for v=1:30 %this sets the outer loop, which ramps rainfall per time step from well below potential ET to well above it

% in small increments

sumCaconc(v)=0;% this and the following lines reset summed and mean variables, so that the run at each rainfall level is

% independent of previous runs.Variables explained where they appear in the

% runs.

sumCamineral1(v)=0;

sumCamineral2(v)=0;

sumCamineral3(v)=0;

sumtotalrockCa(v)=0;

sumproportiontotalCarem(v)=0;

sumweatherCa(v)=0;

sumCaleach(v)=0;

sumrain(v)=0;

meanCaconc(1)=1;

meanCamineral1(1)=0;

meanCamineral2(1)=0;

meanCamineral3(1)=0;

meantotalrockCa(1)=0;

meanproportiontotalCarem(1)=1;

meanrain(v)=0;

waterm(v)=v;% sets rainfall for a given iteration (of the outer loop). Where rainfall does not have a random

%component, this is rainfall; where it does have a random component,

%this is the starting point for calculating rainfall.

% set first elements in arrays, which initializes the initial conditions of

% cations to the same point for each rainfall level.

Camineral1(1)=1000; %pool size of cations in mineral 1

Camineral2(1)=1000; %cations in mineral 2

Camineral3(1)=1000; %cations in mineral 3

totalrockCa(1)=3000; %total cations in primary minerals in the soil

soilwatercarryover(1)=0; % the following terms are still initializations; will be explained when they are used in

% calculations

rain(1)=waterm(v);

soilwaterpostrain(1)=0;

actET(1)=0;

soilwaterpostET(1)=0;

Caconc(1)=0;

Caleach(1)=0;

waterleach(1)=0;

soilwaterpostleach(1)=0;

Caconcpostleach(1)=0;

weatherCa1(1)=0;

Caconcpost1(1)=0;

weatherCa2(1)=0;

Caconcpost2(1)=0;

weatherCa3(1)=0;

Caconcpost3(1)=0;

totalweather(1)=0;

proportiontotalCarem(1)=1;

year(1)=1;

% begin loop in time steps, where calculations take place

for y=begin:stop

x=y-1;

year(y)=y;

ifvariablerain<1 % if there is no random component to rainfall

rain(y)=waterm(v);

else

rain(y)=waterm(v)+(waterm(v)*(10*rand-6.5)); % a random component to rainfall

end

if rain(y)<0 % makes sure rain can't be negative

rain(y)=0.001;

end

if y>changerain % if rainfall changes systematically at a particular point during a run, that is introduced here.

rain(y)=rain(y)*changeamnt;

end

soilwaterpostrain(y)=soilwaterpostleach(x) + rain(y); %takes water left over from last year (after ET and leaching)

%and adds rainfall to give total water

ifsoilwaterpostrain(y)<WHC %if that total is less than water holding capacity

actET(y)=(soilwaterpostrain(y)/WHC)*PotET; %then actual water loss by ET is a fraction of potential ET (accounts

% for water limitation to ET)

else

actET(y)=PotET; % otherwise actual water loss by ET is equal to potential

end

soilwaterpostET(y)=soilwaterpostrain(y)-actET(y); %water left in soil after rain and ET

waterleach(y)=soilwaterpostET(y)-WHC; % water lost by leaching is then water left in soil minus water holding

% capacity; the vision here is that rain adds water in multiple small

% increments during the time step, so ET has a chance to work before

% leaching. If rain came in one big dump, these should be reversed.

ifwaterleach(y)<0

waterleach(y)=0;

end

Caleach(y)=(waterleach(y)/soilwaterpostET(y))*Caconc(x); % this step leaches cations from the soil, Caconc is the

% total soluble cation in the soil, and a proportion of that leaches

% that is equivalent to the proportion of water in the soil that

% leaches.

Caconcpostleach(y)=Caconc(x)-Caleach(y); %soluble cations left in the soil

ifCaconcpostleach(y)<0

Caconcpostleach(y)=0

end

soilwaterpostleach(y)=soilwaterpostET(y)-waterleach(y); % the soil water pool after rain, ET, and leaching

if Camineral1(x)<2 %if mineral 1 is almost gone (here .2% of initial), then it doesn't weather further

weatherCa1(y)=0;

elseifCaconcpostleach(y)>MaxCaconc1 %or if the soluble cation pool after leaching is greater than the solubility of

%the mineral, it doesn't weather

weatherCa1(y)=0;

else

weatherCa1(y)=weathercoeff1*(soilwaterpostleach(y)/WHC)*((MaxCaconc1-Caconcpostleach(y))); % otherwise it

%weathers at a rate proportional to its kinetic coefficient times

%the water content of the soil (as a fraction of WHC) times the

%extent of solution disequilibrium (maximum solubility minus actual

%concentration).

end

Caconcpost1(y)=Caconcpostleach(y)+weatherCa1(y); %cations in solution after mineral 1 has weathered;

%assumes most soluble mineral weathers first.

if Camineral2(x)<2 %weathering of next most soluble mineral is calculated similarly

weatherCa2(y)=0;

elseif Caconcpost1(y)>MaxCaconc2

weatherCa2(y)=0;

else

weatherCa2(y)=weathercoeff2*(soilwaterpostleach(y)/WHC)*((MaxCaconc2-Caconcpost1(y)));

end

Caconcpost2(y)=Caconcpost1(y)+weatherCa2(y); %concentration in solution after both minerals 1 and 2 have

%weathered

if Camineral3(x)<2 %weathering of the least soluble mineral is calculated similarly

weatherCa3(y)=0;

elseif Caconcpost2(y)>MaxCaconc3

weatherCa3(y)=0;

else

weatherCa3(y)=weathercoeff3*(soilwaterpostleach(y)/WHC)*((MaxCaconc3-Caconcpost2(y)));

end

Caconcpost3(y)=Caconcpost2(y)+weatherCa3(y); %quantity of soluble cation in soil after all three minerals

% have weathered

Caconc(y)=Caconcpostleach(y)+weatherCa1(y)+weatherCa2(y)+weatherCa3(y); %soluble cations in solution after weathering

%calculated as the residual after leaching plus weathering release from

%all three minerals.

Camineral1(y)=Camineral1(x)-weatherCa1(y); %reset the pool of primary minerals, and the total of all three minerals

% in the soil

Camineral2(y)=Camineral2(x)-weatherCa2(y);

Camineral3(y)=Camineral3(x)-weatherCa3(y);

totalrockCa(y)=Camineral1(y)+Camineral2(y)+Camineral3(y);

proportiontotalCarem(y)=totalrockCa(y)/3000; %the fraction of the initial primary mineral cations (all three minerals

%that remains in primary form

sumrain(v)=sumrain(v)+rain(y); % calculates the sum of rainfall received, across all the time steps

if y>=startaverage% defines the point where results begin to be summarized; useful for some of these variables

% because they are sensitive to random variations in rainfall.

% Others aren't sensitive in this way (like residual minerals) but

% all are calculated the same

sumCaconc(v)=sumCaconc(v)+Caconc(y); %sum across years at a given rainfall for soluble cations

sumCamineral1(v)=sumCamineral1(v)+Camineral1(y); %sum for each mineral, for weathering, for total primary mineral

%cations, and for the proportion of initial cations remaining in

%primary minerals

sumCamineral2(v)=sumCamineral2(v)+Camineral2(y);

sumCamineral3(v)=sumCamineral3(v)+Camineral3(y);

sumweatherCa(v)=sumweatherCa(v)+weatherCa1(y)+weatherCa2(y)+weatherCa3(y);

sumtotalrockCa(v)=sumtotalrockCa(v)+totalrockCa(y);

sumproportiontotalCarem(v)=sumproportiontotalCarem(v)+proportiontotalCarem(y);

end

end

meanCaconc(v)=sumCaconc(v)/(stop-startaverage); % mean quanity of soluble cations in the soil, for the time increment

%over which final results are determined

meanCamineral1(v)=sumCamineral1(v)/(stop-startaverage); % mean primary mineral cations remaining, and proportions, and

%total weathering, for the time increment.

meanCamineral2(v)=sumCamineral2(v)/(stop-startaverage);

meanCamineral3(v)=sumCamineral3(v)/(stop-startaverage);

meantotalrockCa(v)=sumtotalrockCa(v)/(stop-startaverage);

meanproportiontotalCarem(v)=sumproportiontotalCarem(v)/(stop-startaverage);

meanweather(v)=sumweatherCa(v)/(stop-startaverage);

meanrain(v)=sumrain(v)/(stop-begin); %mean rainfall over the entire run

graph(1,v)=waterm(v); % fill in the output array; this includes the final state of the system at or near the end of the

%time step runs, for each rainfall considered

graph(2,v)=meanrain(v);

graph(3,v)=meanCaconc(v);

graph(4,v)=meanproportiontotalCarem(v);

graph(5,v)=meanweather(v);

end

outgraph=transpose(graph); %transposes the output matrix, and saves it in a file that is readable in excel

savegrad.txtoutgraph-ascii-tabs

figure %calls a figure, here plotting mean rain through the whole run (for each rainfall level) as x axis and the mean

%of soluble cations (which I assume to be in equilibrium with exchangeable

%cations), of the proportion of primary mineral cations remaining as

%primary minerals, and as the ongoing weathering rate at the end of the run

plot (meanrain,meanCaconc,'k', meanrain,meanproportiontotalCarem,'k:', meanrain,meanweather, 'r')

gridoff

axis ([0 32 0 1.1])

ylabel ('Ca concentration and Ca remaining')

xlabel ('rainfall per time step')