\section{Estimating MASP vehicular emissions using VEIN model}

This section presents the application of the most important functions of the VEIN model. These functions obtain an estimate of CO emissions from LDV fleets in MASP for 2015 (for a typical non-holiday week).


\subsection{Traffic data for MASP}	
Hourly traffic is a requirement for this data. This data can only be represented as one hour of data, which can then be extrapolated with the VEIN functions or as a list of hourly traffic data for all the covered hours. The present application includes a \textcolor{red}{4-stage travel demand model} for MASP \cite{CET2014} for a morning rush hour loaded into R as a {\tt SpatialLinesDataFrame}. It includes peak and free flow speeds, along with capacity (maximum amount of vehicles that can circulate in a road per hour) and traffic flow from LDV and HDV. Fig. (\ref{fig:ldv}) shows the traffic simulation of LDV at 08:00-09:00 local time (LT), where urban motorways concentrate the highest amount of vehicles. The total volume of LDV is 24708767 $veh \cdot h^{-1}$ and the number of streets are 34733 with mean of 711 $veh \cdot h^{-1} \cdot street^{-1}$. It is important to note that the VEIN model provides an extraction of the traffic simulation for the western part of São Paulo. The traffic simulation for MASP has a size of 61.6 Mb and the extraction for the western part of São Paulo 3.1 Mb. We provided the extraction and not the whole traffic simulation in VEIN to make it faster. This section provides codes to run VEIN so that the reader can follow them with the data provided in the model. 

	\begin{figure*}
		\centering
		\includegraphics[width=15cm]{figures/ldv.png}
		\caption{Traffic flow simulation for LDV ($\mathrm{veh \cdot h^{-1}}$) at 08:00-09:00 LT for MASP.}
		\label{fig:ldv}
	\end{figure*}

After loading traffic data, the traffic flow was expanded to each hour of the week with the function \textit{temp\_fact}, as shown in the following scripts. It is also necessary to extrapolate hourly vehicle speeds. Therefore, we created the function \textit{netspeed}, which applies the function BPR \citep{bpr} curves, according to the Eq. \ref{eq:bpr}. To use BPR, a data-frame is required with total traffic at all hours and the morning rush parameters capacity, peak speed, free flow speed, length of the road, and with BPR parameters alpha and beta. The argument \textit{scheme} produces a 24 hour speed data-frame, based only on peak and free flow speed with a profile of free flow speeds at early mornings, peak speeds, morning and evening rush hours, and the average at hours in between. If the time-lapse for the emissions estimation is longer than a week, the user could simply replicate the hours until it reaches the desired hours. 

	\begin{verbatim}
data(net)
data(pc_profile)
pcw <- temp_fact(net$ldv+net$hdv, pc_profile)
speedspeed <- netspeed(pcw, net$ps, net$ffs, net$capacity, net$lkm, alpha = 1)
	\end{verbatim}

For example purposes, the resulting speeds can be observed in Fig. \ref{fig:mapspeed}, which shows two different speed maps: one for 08:00 LT (left panel) and the other for 23:00 LT (right panel). This Figure shows that the highest speeds are found in most of the streets further away from the MASP center at both times (08:00 and 23:00 LT). The major difference between the two panels (left and right), is that in the late night, the flow is faster near the center of MASP. This seems reasonable since the vehicular flow tends to diminish during the night. The average speeds also show a pattern related to the type of street as shows Fig. \ref{fig:speed168}. The type of streets comes form São Paulo traffic simulation and they were translated to English and are defined as Motorway: roads with speed limits above 80 $km \cdot h^{-1}$ without physical intersections. Arterial: roads with speed limit of 60 $km \cdot h^{-1}$ with intersection such as traffic lights. Collector: roads with speed limit of 40 $km \cdot h^{-1}$ that collect and distribute traffic between Arterial streets. Local: roads with speed limit of 30 $km \cdot h^{-1}$ that access restricted zones. Figure \ref{fig:speed168} shows that lower speeds are found during the morning (07:00-10:00 LT) and evening (17:00-20:00 LT) rush hours. This is important in terms of air pollution because at lower speeds vehicles emit more pollutants \citep{NtziachristosSamaras2016}. On the contrary, maximum average speeds for each type of road are obtained at night hours and on Sundays at all hours. 

	\begin{figure*}
		\centering
		\includegraphics[width=18cm]{figures/speed1.png}
		\caption{Traffic speeds (colored lines; $\mathrm{km \cdot h^{-1}}$) for LDV fleet at 08:00 LT (left panel) and 23:00 LT (right panel) in MASP.}
		\label{fig:mapspeed}
	\end{figure*}

	\begin{figure}
		\centering
		\includegraphics[width=15cm]{figures/speed2.png}
		\caption{Traffic average speeds for LDV fleet by type of street (colored lines) at 08:00-09:00 LT in MASP.}
		\label{fig:speed168}
	\end{figure}

    After calculating the São Paulo traffic flow average speeds for each hour of the week and each street link, the age distribution of the fleet was obtained by type of vehicle. The \textit{age*} functions (\textit{age\_ldv}, \textit{age\_hdv} and \textit{age\_moto})  distribute the traffic data by the vehicle's age of use. These functions return a data-frame with the number of rows matching the number of streets, columns representing the amount of vehicles by age of use, and a message indicating the average age of the fleet. 
The \textit{age*} functions are related to the Eq. \ref{eq:qvca}, where they split the vehicular flow at street link $Q_i$ by type of vehicle $j$ and the vehicle's age of use $k$. These functions are based on the survival equations presented in the Brazilian Emissions Inventory Report by the \cite{MMA2011} and parameterized for the VEIN model. They allow the use of different coefficients to obtain different age distributions allowing the representation of different realities. Furthermore, the function \textit{my\_age} distributes the traffic from an existing dataset, e.g. yearly vehicle licensing. 
    
    The following code shows three uses of \textit{age*} functions. The first, \textit{my\_age}, uses yearly traffic data from the report of the São Paulo emissions inventory \citep{CETESB2015}, expressed as CETESB\_PC based on vehicle sales. The second, \textit{age\_ldv}, uses default parameters, and the third, \textit{age\_ldv}, uses b = -0.14. 
    
		\begin{verbatim}
CETESB_PC <- c(33491,22340,24818,31808,46458,28574,24856,28972,37818,49050,87923,
133833,138441,142682,171029,151048,115228,98664,126444,101027, 84771,55864,36306,
21079,20138,17439, 7854,2215,656,1262,476,512, 1181, 4991, 3711, 5653, 7039,5839,
4257,3824,3068)
pc1 <- my_age(x = net$ldv, y = CETESB_PC, name = "PC")
pc2 <- age_ldv(x = net$ldv, name = "PC", agemax = 41)
pc3 <- age_ldv(x = net$ldv, name = "PC", b = - 0.14, agemax = 41)
		\end{verbatim}
		
        Figure \ref{fig:Figura1} (a) shows three age distributions, each one sums 24708767 $\mathrm{veh \cdot h^{-1}}$, and each have different average age. It represents São Paulo LDV fleet according to its age of use, with an estimated average age of 11.09 years (red line), 15.53 years (blue line), and 15.17 years (green line ). 
 This Figure shows that in MASP there are more newer vehicles than older ones. 
\textit{age*} functions also include a logical option named \textit{bystreet}, with default value equal to FALSE. When this value is TRUE, \textit{age*} expects that the coefficients \textit{a} and \textit{b} for the \textit{age*} functions are numeric vectors, with length matching the number of streets. This allows different age distributions within the same road network and it is particularly useful for areas %inside a city restrictions.
with less information about the vehicles' age of use. 

\subsection{Emission Factors}

Once we obtain the traffic flow for the desired type of vehicles (in our LDV example), for each hour of the day, for all (desired) days of the week, for each age distribution and for each street link, then we can proceed to the emissions calculation by itself.

The VEIN package includes a database titled \textit{fe2015} with emission factors for PC and light trucks by age of use from the São Paulo official vehicular emissions inventory \citep{CETESB2015}. This inventory was performed using a top-down approach %The method of this inventory is top-down \citep{CETESB2015}
 and the pollutants estimated were CH$_4$, CO, CO$_2$, HC, N$_2$O, NMHC, NO$_X$, and PM. This data includes national and equivalent Euro Emission Standards by year and age. The equivalence among Brazilian \cite{PROCONVE}, \cite{Promot}, \cite{Promo} and Euro \cite{euro} was added into this database in order to choose the corresponding matching vehicle and emissions standard. The equivalence can be seen in Table \ref{tab:euro}.
    
\input{tables/tab_euro.tex}

    \textit{fe2015} emission factors do not include the deterioration effect due to accumulated age of vehicle and it must be included. This is performed with the the deterioration factor function \textit{emis\_det} which has the arguments: pollutant, size of engine, Euro Standard and mileage in km. VEIN includes a Brazilian database of mileage functions named \textit{fkm}, which is a list of functions with each element of the list corresponding to vehicle type. These functions depend on the vehicle's age of use and they originate from the odometer readings of more than $\mathrm{1.6 \cdot 10^6}$ vehicles \citep{BruniBales2013}. 
    
Emission factors for PC, LCV and motorcycles are called with the function \textit{ef\_ldv\_speed}. In the case of trucks and buses they use the function \textit{ef\_hdv\_speed}. %These functions behave as a look-up table. 
    The arguments are filters for an internal database of emission factors which include several parameters such as fuel, Euro Standard, volume of engine and load, among others. These functions also include a multiplication factor with a default value of 1. Exact spelling is required when using the arguments. If the argument names are entered incorrectly, VEIN will not return the emission factor functions. 
    
    The following code shows how to read the emission factors of the VEIN databases $fe2015$, $pc\_profile$ and $fkm$, in order to incorporate the deterioration effect into the \cite{CETESB2015} emission factors. The age of LDV shown in Fig. \ref{fig:Figura1} (a) has a length of 41 years. This means that it needs 41 emission factors, one per each age of use. \textcolor{red}{ (14) We used a 41 years distribution because we want to have a realistic representation of the vehicles in circulation by age of use, despite that oldest vehicles have a small number. This allow us to obtain emissions by age of use of the vehicles, hence, VEIN can be used as a fast tool for scrapping policies and other applications}. It calls the function \textit{emis\_det}, which requires the accumulated mileage, obtained from the list of mileage equations $fkm$. The Fig. \ref{fig:Figura1}(c) shows the emission factors from CETESB with and without deterioration by age of use. We are using deterioration factors from \cite{NtziachristosSamaras2016} that affect only vehicles with a catalytic system. The base year of this emissions estimation is 2015 and the vehicles with catalytic system started in 1992 in Brazil (23 years before 2015). Therefore, the vehicles that entered into the market before 1992 do not include deterioration. The emission factors dataset $fe2015$ includes emission factors for vehicles with only 36 years of use but the vehicular distribution calculated in the last script has 41 years of use. Therefore, we repeated the oldest emission factors to have 41 emission factors. Here we are assuming that the emission factors of vehicles 36 years of use it is the same as the vehicles till 41 years of use. The last line of the following script calculates the deteriorated emission factors of Passenger Cars by age of use. 
    \begin{verbatim}
data(fe2015)
data(pc_profile)
data(fkm)
pckm <- fkm[[1]](1:24)
pckma <- cumsum(pckm)
cod1 <- emis_det(po = "CO", cc = 1000, eu = "III", km = pckma[1:11])
cod2 <- emis_det(po = "CO", cc = 1000, eu = "I", km = pckma[12:24])
co1 <- fe2015[fe2015$Pollutant == "CO", ]
co1[37:41, ] <- co1[36, ]
cod <- c(co1$PC_G[1:24] * c(cod1, cod2), co1$PC_G[25:nrow(co1)])
	\end{verbatim}

Once the deterioration effect was added into the Brazilian emission factors \cite{CETESB2015}, they were scaled to account for speed with the function \textit{ef\_ldv\_scaled}. 
This function is used to multiply emission factors from \textit{ef\_ldv\_speed} with a constant.

The new emission factor (dependent on speed) has the same value as the local emission factor, which is evaluated at the reference speed of the measurement 34.12 $\mathrm{km \cdot h^{-1}}$ for FTP-75. 
    The default speed value is 34.12 $\mathrm{km \cdot h^{-1}}$, but this value must change correspondingly to the speed of the driving conditions. To use this function, it is necessary to scale emission standards of local emission factors with Euro Standards. In the following code, \textit{Euro\_LDV} is a vector indicating Euro Standard by age of use.
	
	\begin{verbatim}
lef <- ef_ldv_scaled(co1, cod, v = "PC", cc = "<=1400", f = "G", p = "CO", 
eu = co1$Euro_LDV)
	\end{verbatim}
	

