% Nfixwater. Last modified 5/11/98. This file is designed to evaluate the effects
% of rainfall and rainfall variation on N inputs (including fixation) and outputs,
% on 15N natural abundance, and on the openness ofN cycling across moisture
gradients.
% It is modified extensively from Nfix15, which the Biogeochemistry Nfix/CO2
paper
% is based on, especially to make sure production and decomposition occur
% in the same simulated year.
% This is the list of constants
CNcrit = 20;
CNf = 33.33;
CNnf = 50;
CPcrit=200;
fixcost = 0.3;
fixlitinput = 0.1;
% Change fixturnover to alter fraction sent to litter pool each year
fixturnover = 0.1;
floatingCNnf = 50;
NFlitinput = 0.1;
% Change NFturnover to alter fraction sent to litter pool each year
NFturnover = 0.1;
Nlossfactor=0; %for don or transformation dependent n losses
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=1;
yrs=6000; %number of years in run
% This is the list to set the value of the first element in every array.
%for shadeswitch=1:2
% shadecontrol=shadeswitch-1;
% for donon=1:2
% Nlossfactor=(donon-1)*.03;
waterm(1)=0;
mean_non_fixer(1)=0;
meanfixer(1)=0;
mean_delta_N15(1)=0;
%for v= 1:13
v=5;
waterm(v)=(2+v)/10; %mean annual water input, with 1=no water limitation
% This is the list to set the value of the first element in every array.
water(1)=0;
water1(1)=0;
water2(1) = 0;
Pweather(1)=1;
Plab(1)=0;
soilP(1)= 0.0001;
Pmin(1)=0;
minP(1)=0;
Ploss(1)=0;
NPPNF(1) = 0;
NPPfix(1) = 0;
soilC(1) = 0;
soilN(1) = 0.0001;
N15pool(1) = soilN(1)*.003663;
N15abundance(1) = 0.003663;
N15input(1) = 0;
N15fix(1) = 0;
N15don(1) = 0;
N15nitrate(1) = 0;
totalN(1) = .0001;
totalNcheck(1) = .0001;
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;
startminN(1)=0;
midminN(1)=0;
finalminN(1)=0;
potuptake(1) = 0;
realup(1)=0;
masslimitNF(1) = 0;
masslimitF(1) = 0;
shade(1)=0;
Ninput(1) = 0;
year(1)=1;
totalNF(1)=0;
totalfix(1)=0;
totalDN15(1)=0;
total(1)=0;
total15(1)=0;
for y = 2:yrs
x=y-1;
year(y)=y;
% Here we deal with year to year variation. waterm + waterm*(rand-.5) gives
an
% even distribution between half and half-again the annual mean (constrained
to 1 if
% over 1); dividing the second term by .5+waterm reduces variance for increasing
% rainfall. water1 can be used to make decomposition and weathering less sensitive
% to ppt than is production; just set water1 to water^.5. water2 can be used
to make
% leaching losses less sensitive to water than production or decomposition,
set
% water2 to water^2.
if y>=4000
water(y)=waterm(v)+(waterm(v)*(rand-.5)); %even distribution of random rain
% for a normal distribution instead, substitute this after waterm* randn/3)/(.5+waterm);
if water(y)<=0 %relevant only for the normally distributed rain case
water(y)=.1;
end
water1(y)=water(y);
water2(y)=water(y);
else
water(y)=waterm(v);
water1(y)=water(y);
water2(y) = water(y);
end
if water(y)>1
water(y)=1;
water1(y) = 1;
water2(y) = 1;
end
Pweather(y)=10*water1(y)/sqrt(y);
if Pweather(y) <= 2*water1(y)
Pweather(y) = 2*water1(y);
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*water(y);
% 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))^NFBiomass(x))/(100+(1+(.00025))^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;
masslimitF(y) = fixBiomass(x)/5 + 500;
% Nadd lets us pulse an N addition in any year
% There is no addition if the year is out of the simulation range
if y >= 6000
if y <= 6200
Nadd(y) = 40;
else
Nadd(y) = 0;
end
else Nadd(y) = 0;
end
% decomp = 5% of soil C
soilC(y) = soilC(x)+NFlit(x)+fixlit(x);
decomp(y) = soilC(y) * 0.05*water1(y);
soilC(y) = soilC(y) - decomp(y);
soilN(y) = soilN(x)+NFlitN(x)+fixlitN(x);
soilCN(y)=soilC(y)/soilN(y);
soilP(y) = soilP(x)+NFlitP(x)+fixlitP(x);
soilCP(y) = soilC(y)/soilP(y);
% Calculate P mineralization
Pmin(y) = soilP(y) - soilC(y)/CPcrit;
if Pmin(y) <= 0
Pmin(y) = 0.0001;
end
soilP(y) = soilP(y) - Pmin(y);
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.000001;
end
soilN(y) = soilN(y) - (1+.5*Nlossfactor+.5*Nlossfactor*water2(y))*Nmin(y);
startminN(y) = finalminN(x) + Nmin(y) + Nadd(y) + Ninput(y);
% N15 of don outputs (assumed to be non-fractionating, and so at 15N abundance
of the system) and
%of non-fixation inputs (assumed to be 5 per mil depleted).
N15don(y) = (.5*Nlossfactor+.5*Nlossfactor*water2(y))*Nmin(y)*N15abundance(x);
N15input(y) = Nadd(y)*.003663+Ninput(y)*(.003663-.005*.003663);
% Calculate potential N uptake for non-fixer
if startminN(y) <= 0
startminN(y) = .00001;
end
if startminN(y) >= 100*water(y)
potuptake(y) = 100*water(y);
else
potuptake(y)=startminN(y);
end
if potuptake(y) < startminN(y)/2
costNF(y) = 0.1;
else
u = potuptake(y)/startminN(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)*(1-costNF(y));
if masslimitNF(y) <= NPPNF(y)
if masslimitNF(y)/CNnf <startminN(y)/2
costNF(y) = 0.1;
else
u = (masslimitNF(y)/CNnf)/startminN(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
midminN(y) = startminN(y)-Nuptake(y);
minP(y) = minP(y) - Pupt(y);
%Calculate NPP of the fixer
%P limited growth
potPNPPfix(y) = minP(y)*CPf;
%growth potential remaining after growth of non-fixers (limited by light or
water
%on site)
potNNPPfix(y) = 5000*water(y)*(1-fixcost)*(4500*water(y)-NPPNF(y))/(4500*water(y));
% This function calculates the shade limited fixer NPP
potshNPPfix(y) = 5000*(1-fixcost)*(1-shade(y));
%choose the minimum of the NPP estimates for fixers
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 off, convert the following remark to code.
NPPfix(y) = 0;
%calculate quantity of N fixed, amount of 15N fixed (assuming no fractionation,
%so zero per mil), and P taken up by fixers.
Nfixation(y) = NPPfix(y)/CNf;
N15fix(y) = Nfixation(y)*.003663;
minP(y) = minP(y) - NPPfix(y)/CPf;
% Calculate litter production and biomass (before any fire occurs)
NFlit(y) = NFturnover*NFBiomass(x);
fixlit(y) = fixturnover*fixBiomass(x);
NFlitN(y) = NFlitinput*NFBiomassN(x);
fixlitN(y) = fixlitinput*fixBiomassN(x);
NFlitP(y) = NFlitinput*NFBiomassP(x);
fixlitP(y) = fixlitinput*fixBiomassP(x);
NFBiomass(y) = NFBiomass(x) - NFlit(y) + NPPNF(y);
fixBiomass(y) = fixBiomass(x) - fixlit(y) + NPPfix(y);
NFBiomassN(y) = NFBiomassN(x) - NFlitN(y) + Nuptake(y);
fixBiomassN(y) = fixBiomassN(x) - fixlitN(y) + Nfixation(y);
NFBiomassP(y) = NFBiomassP(x) - NFlitP(y) + Pupt(y);
fixBiomassP(y) = fixBiomassP(x) -fixlitP(y) + NPPfix(y)/CPf;
floatingCNnf = NFBiomass(y)/NFBiomassN(y);
% This lets us lose some fraction of excess N mineralized
Nloss(y) = midminN(y)*water2(y);
finalminN(y) = midminN(y) - Nloss(y);
% N15 of losses by denitrification or nitrate leaching, assumed to fractionate
10 per
% mil relative to total system N.
N15nitrate(y) = Nloss(y)*(N15abundance(x)-.01*N15abundance(x));
%summation of N balances - of inputs versus outputs for total N and N-15, and
of the
% within system pools for total N.
totalNcheck(y) = totalNcheck(x) + Ninput(y) + Nfixation(y)-Nloss(y)-...
(.5*Nlossfactor+.5*Nlossfactor*water2(y))*Nmin(y);
totalN(y) = soilN(y) + NFBiomassN(y) + fixBiomassN(y) + finalminN(y)...
+ NFlitN(y) + fixlitN(y);
N15pool(y) = N15pool(x) + N15input(y) + N15fix(y) - N15don(y) - N15nitrate(y);
N15abundance(y) = N15pool(y)/totalN(y);
DN15(y) = 1000*(N15abundance(y)-.003663)/.003663;
% P losses and within-system fractions
Ploss(y) = .1*water2(y)*minP(y);
Plab(y) = Plab(x) + .7*minP(y);
minP(y) = .2*minP(y) + .1*minP(y)*(1-water2(y));
%get a summation of fixer, non-fixer, N15 for the last 1000 years of simulation
for biomass,
%300 yrs for N15.
% if y>(yrs-2000)
% for y1=2:2000
% x1=y1-1;
% totalNF(y1)=totalNF(x1)+NFBiomass(y);
% totalfix(y1)=totalfix(x1)+fixBiomass(y);
% total(y1)=total(x1)+1;
% end
end
% if y>(yrs-300)
% for y2=2:300
% x2=y2-1;
% totalDN15(y2)=totalDN15(x2)+DN15(y);
% total15(y2)=total15(x2)+1;
% end
% end
%end
% print mean biomass of fixer and non-fixer, and system delta 15N, in last
% 300 years of simulation.
% mean_non_fixer(v)=totalNF(y1)/total(y1);
% meanfixer(v)=totalfix(y1)/total(y1);
% mean_delta_N15(v)=totalDN15(y2)/total15(y2)
% plot water input
%figure
%plot (year,water,'b')
%grid on
%axis ([1850 1900 0 1])
%ylabel ('water input')
%xlabel ('year')
% Plot biomass on a linear scale
figure
plot(year,NFBiomass, 'b', year, fixBiomass, 'r--')
grid on
axis([5000 6000 15000 40000])
ylabel ('Biomass')
xlabel('Year')
print biomass1 -dpsc
%plot npp
%figure
%plot (year,NPPNF,'b',year,NPPfix,'r--')
%grid on
%axis ([0 2000 0 5000])
%ylabel ('NPP')
%xlabel ('year')
% 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 Nmin
%figure
%plot(year,Nmin,'b')
%grid on
%axis ([1850 1900 0 100])
%ylabel ('Nmin')
%xlabel ('year')
%plot n loss
%figure
%plot (year, Nloss,'b')
%grid on
%axis ([5000 6000 0 50])
%ylabel ('nloss')
%xlabel ('year')
%plot n15
%figure
%plot (year, DN15, 'b')
%grid on
%axis ([0 2000 -6 12])
%ylabel ('delta 15n')
%xlabel ('year')
%plot totaln
%figure
%plot (year, totalN, 'b', year, totalNcheck, 'r--')
%grid on
%axis ([0 3000 0 5700])
%ylabel ('totaln')
%xlabel ('year')
%end
%plot biomass versus rain
%figure
%plot (waterm, mean_non_fixer, 'b', waterm, meanfixer, 'r--')
%grid on
%axis ([.3 1.5 0 50000])
%ylabel ('biomass')
%xlabel ('mean ppt')
%plot N15 versus rain
%figure
%plot (waterm, mean_delta_N15, 'b')
%grid on
%axis ([.3 1.5 -6 12])
%ylabel ('delta N15')
%xlabel ('mean ppt')
%end
%end