Nitrogen Fixation Model
% This is the version modified by Peter Vitousek, changed from the
% earlier version also called Nfix15.
% Nfix 15 -- 28 aug 98: This is modified from Nfix14 to
% correct a couple of problems with the handling of grazing and fire.
% It is also modified to provide easier parameter adjustment.
% Nfix14 -- 4 October 1997
% this is modified from Nfix4 to include N addition (as a pulse)
% and to allow N loss
% this is mdified from Nfix5 to add soil P
% (and modified from Nfix13) to correctly account for P uptake
%modified from Nfix10 to allow independent handling of N, P, and C
%modified from Nfix12 to adjust mass limit, to avoid double correction
% in cost accounting, and to prevent fixers from entering through
% an established canopy of non fixers
% For Nfix13, I also raised the fixcost coefficient to 0.3
% this is modified from Nfix13 to deposit N and P from grazing in litter
% it is modified from Nfix13 to lose N but deposit P following a fire.
% It is also modified from Nfix13 to make grazers prefer N fixers
% over non fixers, and to make the shade function logistic.
% This is the list of constants
CNcrit = 20;
CNf = 33.33;
CNnf = 50;
%CO2 = 1 for current conditions and some multiple of 1 for the future.
% A multiple of X inplies an NPP stimulation of X for the fixer.
% The multiple is 1/2 X for the non-fixer, though only masslimit
% and the c:element ratio of nonfixers are adjusted.
% The CO2 effect is proportional on maximum NPP, masslimit, shade, fixation
% cost,
% and on the non-fixer NPP necessary to suppress the fixer NPP.
% as of 10/17, the CO2 effect is set in the loop at the begining of the
% execution.
CO2 = 1;
CO3 = (CO2 - 1)/2 + 1;
CPcrit=200;
fireloss = 1;
% Make fireon = 1 for a fire every 20 years after year 300.
% Make fireon = 0 for no fires.
fireon=0;
fixcost = 0.3;
fixgrazemultiple = 5;
fixlitinput = 0.1;
% Change fixturnover to alter fraction sent to litter pool each year
fixturnover = 0.1;
floatingCNnf = 50;
% Alter grazefrac to adjust grazing
grazefrac=0.0;
NFlitinput = 0.1;
% Change NFturnover to alter fraction sent to litter pool each year
NFturnover = 0.1;
Nlossfactor=0.05;
NPnf=10;
NPf=10;
CPnf = CNnf*NPnf;
CPf=CNf*NPf;
% Make Padd to 100 to eliminate P limitation.
% Under normal conditions, Padd can be 0.
Padd=100;
Potprod = 5000;
% For shadecontrol, 1 is with the shade function on.
% 0 is with the shade function off.
shadecontrol = 0;
target = 200;
water=1;
% This is the list to set the value of the first element in every array.
Pweather(1)=1;
Plab(1)=0;
soilP(1)= 0.001;
Pmin(1)=0;
minP(1)=0;
Ploss(1)=0;
NPPNF(1) = 0;
NPPfix(1) = 0;
soilC(1) = 0;
soilN(1) = 0.001;
decomp(1) = 0;
Nmin(1) = 0;
Nloss(1) = 0;
NFlit(1) = 0;
fixlit(1) = 0;
NFlitN(1) = 0;
fixlitN(1) = 0;
NFlitP(1) = 0;
fixlitP(1) = 0;
NFBiomass(1) = 0;
fixBiomass(1) = 0;
NFBiomassN(1) = 0;
fixBiomassN(1) = 0;
NFBiomassP(1) = 0;
fixBiomassP(1) = 0;
Nfixation(1)=0;
Nloss(1) = 0;
Nadd(1) = 0;
minN(1)=0;
potuptake(1) = 1;
realup(1)=0;
masslimitNF(1) = 0;
masslimitF(1) = 0;
shade(1)=1;
NFgraze(1) = 0;
Fgraze(1) = 0;
fire(1) = 0;
Ninput(1) = 2;
burnedNF(1)=0;
burnedF(1)=0;
fixation(1) = 1;
for y = 2:5000
% This loop allows one to increase CO2 at some point during the simulation.
% Rem it out if you want to set the CO2 effect from the parameter list.
if y >= 10000
CO2 = 1.5;
CO3 = (CO2 - 1)/2 + 1;
else
CO2 = 1;
CO3 = (CO2 - 1)/2 + 1;
end
x=y-1;
year(y)=y;
Pweather(y)=10/sqrt(y);
if Pweather(y) <= 2
Pweather(y) = 2;
end
%A small fixed rate of N addition makes it possible for the non-fixer
% to eventually establish, even in the absence of the fixer
Ninput(y) = 2;
% The function shade is intended to make it difficult for the
% fixer to invade when non-fixer biomass is high
% This is different from the NPP function which says that the
% sum of fixer and non fixer NPP is fixed
% This gives a logistic type function that saturates at
%biomass = 50,000 but reaches 90% at biomass = 40,000
shade(y)=((1+(.00025/CO2))^NFBiomass(x))/(100+(1+(.00025/CO2))^NFBiomass(x));
% The idea of the mass limit is that NPP can't jump to the
% maximum in a single year
% With biomass/5 + 500 this can limit only when biomass
% is less than 22500
masslimitNF(y) = NFBiomass(x)/5 + 500*CO3;
masslimitF(y) = fixBiomass(x)/5 + 500*CO2;
%The function fire resets biomass to zero and results of the loss of
%some of the N but none P in the biomass. This fraction of N loss
%is controlled by the parameter fireloss.
if fireon >= 1
if y >= 300
roundyear = round(y/20);
isit = y/20 - roundyear;
if isit ==0
fire(y) = 1;
else
fire(y) = 0;
end
else
fire(y) = 0;
end
else fire(y) = 0;
end
% litter = turnover (typically 10%) of last year biomass,
% plus all of the N and P consumed by grazing,
% plus all of the P lost to fire, plus the possibility of
% some of the N released in a fire.
NFlit(y) = NFturnover*NFBiomass(x);
fixlit(y) = fixturnover*fixBiomass(x);
NFlitN(y) = NFlitinput*NFBiomassN(x) +...
burnedNF(x)*(1-fireloss)/floatingCNnf;
fixlitN(y) = fixlitinput*fixBiomassN(x) + burnedF(x)*(1-fireloss)/CNf;
NFlitP(y) = NFlitinput*NFBiomassP(x) + burnedNF(x)/CPnf;
fixlitP(y) = fixlitinput*fixBiomassP(x) + burnedF(x)/CPf;
% Nadd lets us pulse an N addition in any year
% There is no addition if the year is out of the simulation range
% Convert pulse to zero to prevent it
if y >= 300
if y <= 320
Nadd(y) = 0;
else
Nadd(y) = 0;
end
else Nadd(y) = 0;
end
% decomp = 5% of soil C
decomp(y) = soilC(x) * 0.05;
soilC(y) = soilC(x)+NFlit(x)+fixlit(x)-decomp(x) + 0.5*(NFgraze(x)...
+ Fgraze(x));
soilN(y) = soilN(x)+NFlitN(x)+fixlitN(x)-(1+Nlossfactor)*Nmin(x);
soilCN(y)=soilC(y)/soilN(y);
soilP(y) = soilP(x)+NFlitP(x)+fixlitP(x)-Pmin(x);
soilCP(y) = soilC(y)/soilP(y);
% Calculate P mineralization
Pmin(y) = soilP(y) - soilC(y)/CPcrit;
if Pmin(y) <= 0
Pmin(y) = 0.001;
end
minP(y) = minP(x)+Pmin(y)+Pweather(y)+Padd+0.1*Plab(x);
% Calculate N mineralization
% N mineralized only when N available in excess of
% needs to store soil C at C:N = CNcrit
Nmin(y) = soilN(y) - soilC(y)/CNcrit;
if Nmin(y) <= 0
Nmin(y) = 0.001;
end
minN(y) = minN(x) +Nmin(y) - Nloss(x) + Nadd(x) + Ninput(x);
% Calculate potential N uptake for non-fixer
if minN(y) <= 0
minN(y) = .001;
end
if minN(y) >= 100
potuptake(y) = 100;
else
potuptake(y)=minN(y);
end
if potuptake(y) < minN(y)/2
costNF(y) = 0.1;
else
u = potuptake(y)/minN(y)-0.5;
costNF(y) = (u*(u*.05+0.1))/(.5+u)+(.5*.1)/(.5+u);
end
% Calculate potential P uptake for non-fixer
potPupt(y) = (1/NPnf)*potuptake(y);
if potPupt(y)<=minP(y)
Pupt(y)=potPupt(y);
else
Pupt(y) = minP(y);
end
Nuptake(y)=Pupt(y)*NPnf;
%Calculate NPP of the non fixer
NPPNF(y) = CNnf*Nuptake(y)*CO3*(1-costNF(y));
if masslimitNF(y) <= NPPNF(y)
if masslimitNF(y)/CNnf < minN(y)/2
costNF(y) = 0.1;
else
u = (masslimitNF(y)/CNnf)/minN(y)-0.5;
costNF(y) = (u*(u*.05+0.1))/(.5+u)+(.5*.1)/(.5+u);
end
NPPNF(y) = masslimitNF(y)*(1-costNF(y));
if NPPNF(y) <= 0
NPPNF(y) = 0;
end
Pupt(y) = NPPNF(y)/CPnf;
Nuptake(y) = NPPNF(y)/(CNnf);
end
% After calculating Non fixer NPP, we need to decrement the
% pools of mineral N and P
minN(y) = minN(y)-Nuptake(y);
minP(y) = minP(y) - Pupt(y);
%Calculate NPP of the fixer
potPNPPfix(y) = minP(y)*CPf;
potNNPPfix(y) = 5000*...
CO2*(1-fixcost/CO2)*(4500*CO2-NPPNF(y))/(4500*CO2);
% This function calculates the shade limited fixer NPP
potshNPPfix(y) = 5000*CO2*(1-fixcost/CO2)*(1-shade(y));
if potPNPPfix(y) >= potNNPPfix(y)
NPPfix(y) = potNNPPfix(y);
else
NPPfix(y) = potPNPPfix(y);
end
if NPPfix(y) >= masslimitF(y)
NPPfix(y) = masslimitF(y);
end
if shadecontrol >= 1
if NPPfix(y) >= potshNPPfix(y)
NPPfix(y) = potshNPPfix(y);
end
end
%To turn fixerNPP on, convert the following line to a remark.
% NPPfix(y) = 0;
Nfixation(y) = NPPfix(y)/33.33;
minP(y) = minP(y) - NPPfix(y)/CPf;
% Add grazing
% Biomass is last year*loss to turnover & grazing + current year
% NPP - grazing
if y >= 1
NFgraze(y) = grazefrac*NFBiomass(x);
Fgraze(y) = fixgrazemultiple*grazefrac*fixBiomass(x);
NFlitinput = NFturnover + grazefrac;
fixlitinput = fixturnover + grazefrac*fixgrazemultiple;
else
NFgraze(y) = 0;
Fgraze(y) = 0;
NFlitinput = NFturnover;
fixlitinput = fixturnover;
end
% Calculate Biomass (before any fire occurs)
NFBiomass(y) = (1-NFlitinput)*NFBiomass(x) + NPPNF(y);
fixBiomass(y) = (1-fixlitinput)*fixBiomass(x) + NPPfix(y);
NFBiomassN(y) = (1-NFlitinput)*NFBiomassN(x) + Nuptake(y);
fixBiomassN(y) = (1-fixlitinput)*fixBiomassN(x) + NPPfix(y)/CNf;
NFBiomassP(y) = (1-NFlitinput)*NFBiomassP(x) + Pupt(y);
fixBiomassP(y) = (1-fixlitinput)*fixBiomassP(x) + NPPfix(y)/CPf;
floatingCNnf = NFBiomass(y)/NFBiomassN(y);
% Add fire (if a fire occurs in this year)
if fire(y) == 1
burnedNF(y) = NFBiomass(y);
burnedF(y) = fixBiomass(y);
NFBiomass(y) = 0;
fixBiomass(y) = 0;
NFBiomassN(y) = 0;
fixBiomassN(y) = 0;
NFBiomassP(y) = 0;
fixBiomassP(y) = 0;
else
burnedNF(y) = 0;
burnedF(y) = 0;
end
% This lets us lose some fraction of excess N mineralized
Nloss(y) = minN(y)*water;
Ploss(y) = .1*water*minP(y);
Plab(y) = Plab(x) + .7*minP(y);
minP(y) = .2*minP(y) + .1*minP(y)*(1-water);
end
%Plot NPP
figure
plot(year,NPPNF, 'b', year,NPPfix, 'r--')
grid on
axis([0 500 0 5000])
ylabel('NPP')
xlabel('Year')
print NPP -dpsc
% Plot biomass on a linear scale
figure
plot(year,NFBiomass, 'b', year, fixBiomass, 'r--')
grid on
axis([0 5000 0 50000])
ylabel ('Biomass')
xlabel('Year')
print biomass -dpsc
% Plot biomass on a log scale
figure
semilogy(year,NFBiomass, 'b', year, fixBiomass, 'r--')
%grid on
%axis([0 500 0 50000])
ylabel ('Biomass')
xlabel('Year')
% Plot N and P mineralization
figure
plot(year, Nmin, 'b', year, Pmin, 'r--')
grid on
axis([0 500 0 120])
ylabel('mineralization')
xlabel('Year')
print NPmineralization -dpsc
% Plot N cost to non fixer
figure
plot(year,costNF, 'r')
grid on
axis([0 500 0 0.15])
ylabel('N cost to non-fixer')
xlabel('Year')
% Plot N fixation
figure
plot(year, Nfixation, 'r')
grid on
axis([0 500 0 120])
ylabel('mineralization')
xlabel('Year')
print Nfixation -dpsc