Merge branch 'master'
[wrf-fire-matlab.git] / chem / FEMIS.m
blobbf09e2880fb786a60f2c2b6daab3f7613363a897
1 function [E_CO2, E_CO, E_CH4, E_PM25] = FEMIS(fgi, burn_rate, U_F, V_F, RH4, mode)
2 %FEMIS function computing emission rates of CO2, CO, CH4 and particlate matter
3 %based on:
4 % fgi - fuel load [kg/m2]
5 % burn_rate [kg/s] - mass of the fuel burnt in unit time.
6 % U_F, V_F [m/s]  mid-flame wind speed components (used for computation of the WS_F)
7 % RH4 [%] -relative humidity 4 hours before (?)
8 % mode - 0 - default, computes consumption rates according to FEMIS
9 %           allowing for the total consumption greater than the combustion rate
10 %      - 1 - limits maximum consumption rate to the actual burn rate
11 %      - 2 - computes emissions straight from burnt rate ignoring smoldering adjustment 
13 %all the constants and equations come from the following US Forest Service
14 %technical report: Fire Emission Production Simulator (FEPS) User's Guide
15 %Version 1.0, Gary K. Anderson, David V. Sandberg and Robert A. Norheim
16 %January 2004
17 %available from: www.fs.fed.us/pnw/fera/feps/FEPS_users_guide.pdf
19 %Adam Kochanski 05.26.2012 
21 %CONSTANTS
23 %emission factor constants according to FEPS user's manual p. 63. [g/kg]
24 Kef_CO2=0;
25 Kef_CO=961;
26 Kef_CH4=42.7;
27 Kef_PM25=67.4;
29 %emission factor coefficients according to FEPS user's manual p. 63. [g/kg]
30 EF_CO2=1833;
31 EF_CO=984;
32 EF_CH4=43.2;
33 EF_PM25=66.8;
35 %flam involvment sensitivity coeffcient (eq 5, p. 53)
36 k_agi =1.0;
38 % combustion effciency factors according to FEPS user's manual p. 62 [-]:
39 k_f=0.9;     %flame phase (f)
40 k_sts=0.76;  %short-term smoldering (sts)
41 k_lts=0.76;  %long-term smoldering (lts)
43 %flaming phase consumption coeffcients from to FEPS user's manual p. 53 [-]
44 k_cag=0.5;
45 k_cbg=0.2; % not used since Cbg=0
47 % Duff bulk density [tons/acre-inch] p. 54
48 B_f=20;
49 % STS depth sensitivity to consumption [tons/acre-inch] p.55
50 B_sts=12;
52 %flaming phase resindence time coeficients  (eq. 8, p. 54)
53 k_tflam1= 4/3;
54 k_tflam2= 8.0;
56 %short-term moldering residence time coeffcients (eq.13 p. 55)
57 k_edr1=8/3;
58 k_edr2=0.5;
60 %flame phase diffusivity dependency
61 N_tflam=0.5;
63 %short-term smoldering diffusivity dependency
64 N_edr=0.5;
66 %consumption threshold for flame involvment [tons/acre], 
67 % 1 t/acre = 0.224 kg/m2
68 Cti=10;          % [tons/acre]
69 Cti=Cti * 0.224; % [kg/m2] 
71 %wind speed benchmark for smoldering 1 mph = 0.44704 m/s
72 u_b= 3; %[mph]
73 u_b= u_b * 0.44704; % [m/s]
75 %benchmark relative humidity deficit [%]
76 RHb=60.0; 
78 %Molecular weights [g/mol]
79 MW_CO2=44.0; 
80 MW_CO=28.0;
81 MW_CH4=16.0;
83 %COMPUTATIONS
84     
85     % by default computation from FEPS allowing greater and smaller consumption than the burn_rate
86     
87     %==========================================================================
88     % computing flaming phase consumption, assuming no below ground consumption
89     % Cbg=0; and the above ground consumption (canopy, shrub, grass, wood
90     % litter) Cag equal to fuel combustion rate. eq. 6 p. 53
92     % flaming phase involvment
93     % flaming phase involvment
94     Cag=fgi; % we assume above the ground consuption equal to the fuel load [kg/m2]
95     %Cbg=0;   % we assume no below ground consumption
97     Inv_f=100*(1-k_agi*exp(-Cag/Cti))
99     %==========================================================================
100     % computing short-term smoldering phase consumption [kg/m2]
101     Inv_sts=Inv_f
103     % computing flaming phase consumption rate [kg/s]
104     CR_f= burn_rate * Inv_f/100 * k_cag 
106     % computing smoldering adjustment from eq. 26, p 59
107     WS_F = sqrt (U_F^2 + V_F^2); % wind speed at midflame height [m/s]
108     %U_F midflame X wind speed component
109     %V_F midflame Y wind speed component
110     SA=sqrt(WS_F/u_b) * ((100-RH4)/RHb) 
112     %computing sort-term smoldering consumption rate [kg/s]
113     CR_sts= SA* burn_rate * Inv_sts/100 * k_cag
115     % we assume no long term smoldering (>2h) since in most cases withing
116     % 2h all fuel within a grid cell will burn, and we have no data about the
117     % duff layer
118     CR_lts = 0;
120     %computing total consumption rate according for FEPS (eq. 30, p.61):
121     CR_total=CR_f + CR_sts + CR_lts; %[kg/s]
123     if (mode==1)  % limit the consumption rate to the burn rate
124         if (CR_total>burn_rate) 
125             CR_total=burn_rate;
126         end
127     elseif (mode==2) % compute the consumption rate straight from burn rate
128             CR_total=burn_rate;
129     end
130 CR_total
132 %computing combustion efficiency:
133 CEff=(k_f*CR_f + k_sts*CR_sts + k_lts*CR_lts)/(CR_f + CR_sts + CR_lts + 1)
135 ER_CO2= (Kef_CO2 - EF_CO2*CEff)*CR_total/1000 %[kg/s]
136 ER_CO=(Kef_CO - EF_CO*CEff)*CR_total/1000 %[kg/s]
137 ER_CH4=(Kef_CH4 - EF_CH4*CEff)*CR_total/1000 %[kg/s]
138 ER_PM25=(Kef_PM25 - EF_PM25*CEff)*CR_total/1000 %[kg/s]
140 E_C02 = ER_CO2/(MW_CO2*3600) % CO2 emissions in mol/hr
141 E_C0  = ER_CO /(MW_CO *3600) % CO emissions in mol/hr
142 E_CH4 = ER_CH4/(MW_CH4*3600) % emissions in mol/hr
143 E_PM25 = ER_PM25*1000 % PM25 emissions in ug/s