Derivation of the Euler–Lagrange System for Weak-Constraint 4D-Var Data Assimilation of Simulated Wildfire Smoke PM2.5 into an Eulerian Transport Model

Derivation of the Euler–Lagrange System for Weak-Constraint 4D-Var Data Assimilation of Simulated Wildfire Smoke PM2.5into an Eulerian Transport Model By Sandra Rebecca Babyale ([email protected]w) African Institute for Mathematical Sciences (AIMS), Rwanda Co-supervised by: Professors Jodi Mead and Donna Calhoun Boise State University, USA June 2021 Submitted in partial fulfillment of the requirements of a Master of Science in Mathematical Sciences at AIMS Rwanda Abstract The need to quantify smoke concentration over a large area is vital in forecasting air quality. The pollutant of principle concern in smoke is fine particulate matter with diameters 2.5 micrometers and smaller (PM2.5) because it can enter lungs and blood streams causing respiratory complications and other diseases and it also absorbs radiations causing warming of the atmosphere thus climate change. Complex mesoscale models are used to accurately explain the dynamics and reactions of PM2.5, but these models contain unknown inputs that can be identified using data assimilation. In this study, the weak constraint four-dimensional variational data assimilation (4DVAR) method one of the most comprehensive variational data assimilation techniques is used to assimilate PM2.5data into the one dimensional Eulerian model representing the dynamics and reactions of PM2.5. Here we minimize an approximate cost function that gives the Euler Lagrange systems which are then solved for optimal estimates of PM2.5concentration and other model parameters. i Contents Abstract i 1 Introduction 1 1.1 Background ....................................... 1 1.2 Problem statement ................................... 2 1.3 Aim and objectives ................................... 2 2 Literature Review 3 2.1 Overview of long-range regional transport of smoke ................. 3 2.2 Overview of weak constraint 4D-Var DA technique ................. 3 3 Weak Constraint 4DVAR 6 3.1 The model ........................................ 6 3.2 Data ........................................... 6 3.3 The cost function .................................... 7 3.4 Derivation of Euler Lagrange systems separately ................... 8 3.5 Derivation of the Euler Lagrange system simultaneously .............. 16 4 Pseudo Code for the E-L Systems 20 4.1 Sample pseudo code for ˆq................................ 20 4.2 Sample pseudo code for ˆ Q............................... 20 4.3 Sample pseudo code for ˆ K............................... 22 4.4 Sample pseudo code for ˆc................................ 23 5 Conclusion and Recommendations 24 References 27 ii 1. Introduction 1.1 Background Over the last decade, climate change monitoring has become a fundamental activity due to the different effects of climate change such as global warming, increase in sea level and ocean temperatures, and intense droughts that threaten the growth of crops and so on. This has led to the development of studies of the different factors leading to these changes. Atmospheric particulate is one of the major causes of climate change. These particulate emissions to the atmosphere are generated by both natural processes such as volcanic eruptions, and wildfires, and human activities such as quarrying, burning of fossil fuels, agricultural practices,and many others (Benedetti et al.,2009). According to Ito and Penner (2004), biomass burning, both natural for example wildfires and anthropogenic such as burning for agricultural practices, fuelwood for cooking and heating, to mention, but a few, is a major source of global atmospheric gaseous and particulate emissions (smoke), accounting for approximately 34–38 % and 40% of global carbonaceous aerosol (suspension of fine particles in the air) and black carbon loadings, respectively. Because of their significant radiative effects and climate consequences, aerosol emissions from biomass burning are of particular concern around the world. Smoke particles generated by biomass combustion have a direct effect on air quality, weather, and climate. PM2.5, or small particulate matter with a diameter of 2.5 micrometers or less, is a major concern in smoke because it can reach the lungs and bloodstream, causing health issues like respiratory complications and heart disease (Levine, 1998). Smoke transport, once considered solely a local issue, is now recognized as a complex issue influenced by regional, hemispheric, and even global factors. Despite the fact that domestic sources are the primary causes of much of a given country’s air quality issues, countries are both sources and receivers of contaminants transported over long distances. Pollutants migrate not only from country to country but also among continents (Council et al.,2010). These contaminants cause public health problems, reduced visibility, cause damage to agricultural and natural plants, decreased domestic and wild animal viability, and damage to infrastructure (Steffen et al., 2016). Long-range transport of smoke from foreign sources is thus receiving more attention in both the scientific literature and the popular press. Organizations in charge of ensuring that air quality standards are met are increasingly concerned, making the stakeholders have various advanced viewpoints on how localities should address the impact of international pollutant transport in air quality planning (Bertschi and Jaffe,2005;Council et al.,2010). The need to quantify the concentration of the distant smoke pollutants necessitates the use of complex mesoscale models to accurately explain the processes and dynamics and to adequately integrate chemistry and emissions. For decades, the two components of the related systems, atmospheric physics/dynamics with the largest sub-class of numerical weather prediction models (NWP) and atmospheric chemistry models known as chemistry transport models (CTM), have been developing separately 1 Section 1.2. Problem statement Page 2 for historical and practical purposes . However, the basic features and consistency of the evaluation and forecasting of both the condition and composition of the atmosphere calls for Data Assimilation, a process by which observations of a system are incorporated into the model state of a numerical model of that system (Ghil et al.,1981). Assimilation of observations into numerical models has become a crucial modeling aspect that has led to the development of numerous schemes. These Data Assimilation schemes have developed into extremely complex systems over time, such as the four-dimensional variational method (4DVAR) at the European Center for Medium-Range Weather Forecasts (ECMWF), which is a special case of a generalized least squares problem (Talagrand,1997). This scheme uses a linearized forecast model to ensure that the observations are given dynamically as well as statistically likely responses in the analysis and can handle a wide range of meteorological data. One of the two approaches to 4DVAR is the weak constraint approach which was used for this study. Weak constraint 4DVAR assumes errors in the model’s external forcing, its initial and boundary conditions and hence making the dynamical problem nonlinear even though the model is linear. This approach results into a coupled E-L system of equations hence employing the representer method (Eknes and Evensen,1997). The representer method is used because of its ability to completely decouple the Euler Lagrange system of equations and the solution is obtained in a finite-dimensional space spanned by the representers, one for each measurement (Van Leeuwen and Evensen,1996). 1.2 Problem statement We plan to simulate the long-range regional transport of smoke in one dimension by estimating the concentrations of fine particulate matter with diameters 2.5 micrometers and smaller (PM2.5) in the near surface with a one dimensional transport model on an Eulerian grid given by; qt(x, t) + ∇ · u(x, t)q(x, t)−∇·(K· ∇q(x, t)) + cq(x, t) = Q(x, t),(1.2.1) where q(x, t)is the concentration of PM2.5in (µgm−3), q0(x,t) is the initial concentration of PM2.5 in (µ gm−3), Q(x, t)is the source term which denotes the rate of mass influx per unit volume per unit time into the smoke plume (kgm−3s−1), u(x, t)represents dispersion by ambient wind in (ms−1), Kis a one dimensional diffusivity constant representing turbulent diffusion in (m2s−1) and cis the decay coefficient of the fall out due to rain or temperature in (Bq). 1.3 Aim and objectives The main aim of this project is to simulate PM2.5data into the Eulerian transport model using weak constraint 4DVAR with representers and estimate the errors. The specific objectives are: to derive Euler Lagrange systems of equations when optimizing for q(x, t),Q(x, t),u(x, t),K, and cand decoupling these systems. 2. Literature Review 2.1 Overview of long-range regional transport of smoke Wang et al. (2006) clearly express the role smoke plays in human health and climate change stating that it [smoke] affects the air quality and climate due to an increase in the concentrations of particulate matter and greenhouse gases such as carbondioxide in the atmosphere. Smoke which was once a purely local issue was found by Locatelli et al. (2015) to be an international pollutant of great concern by the air quality management stakeholders who came up with ways on how governments should mitigate the transport of these pollutants to and from their countries. For these to be successful, a quantitative understanding of the magnitudes and dynamics of the flows was needed. Similar research was done by Vijayakumar et al. (2016) who also gave the importance of monitoring the long-range transport of smoke and that the key issue in this is predicting the magnitude and location of smoke effects. Based on these research findings came my zeal to study smoke transport, and more specifically focusing on PM2.5concentrations. However, the underlying objective is how this can be done! Additional findings by Locatelli et al. (2015) suggested the need to quantify the concentration of distant PM2.5using complex mesoscale models to accurately explain the processes and dynamics, as well as to adequately integrate chemistry and emissions. As a result, many experts and academics in different regions have conducted extensive research on PM2.5concentration prediction using both general circulation models (GCMs)with the biggest subset being NWP models that simulate atmospheric dynamics and CTMs which simulate atmospheric chemistry, for practical purposes. Many CTMs have been employed in estimating both ground and atmospheric concentrations of PM2.5, for example, Zhou et al. (2018) used the WRF-Chem to reproduce PM2.5 concentrations although it overpredicted wind velocity and concentrations in certain areas, the GEOS-Chem global chemical transport model was developed by Van Donkelaar et al. (2010) to estimate surface PM2.5concentrations. Besides the CTMs, most studies on long-range simulation of PM2.5have embarked on using NWP models which are governed by a set of time-dependent partial differential equations corresponding to certain conservation laws of physics for example conservation of mass, momentum, transport equations, and on as presented by Ghil et al. (1981). However, their research still suggested that for better estimation of PM2.5concentrations, the models could be used alongside putting chemical reactions into consideration. In their conclusive remarks, they suggested that despite the recent advancements in all these models, they still exhibit errors that need to be reduced and DA was proposed since it provides the means to produce suitable state estimates as it combines models and actual observations. 2.2 Overview of weak constraint 4D-Var DA technique Data assimilation as presented by Xiao and Friedrichs (2014), is an important technique that has been widely used in atmospheric research as it evolved from the area’s rapid theoretical expansion and the desire to translate it into practical applications. They referred to it as a technique 3 Section 2.2. Overview of weak constraint 4D-Var DA technique Page 4 that combines observation of atmospheric conditions with the system’s underlying dynamical principles to produce the best estimate of the system’s state where the model’s uncertainties are reduced than that which could be obtained when using the data and model separately. Ide et al. (1997), argued that the establishment of a global atmospheric data network that could provide the necessary measurements allowed the use of data assimilation techniques since they involve the systematic use of these measurements to constraint a mathematical model. Data assimilation (DA) therefore, provided a tool for using these findings to boost the operational models’ forecasting abilities. There are numerous DA techniques available, each with its own set of limitations. According to Wang et al. (2000), the development of DA methodology has mainly experienced three stages with the first being the simple analysis methods which were the earliest bases of data assimilation and were mostly used in the 1950s before the era of computers. Statistical considerations were introduced into the atmospheric data assimilation in the 1960s and 1970s and some forms of optimum interpolation were used to assimilate observations into forecast models based on these considerations. The last stage, which is being used to date, comprises sequential and variational DA techniques as grouped by Talagrand (1997). Across several studies, Evensen (2002) stated that sequential DA methods which do not require additional features in the modeling process have proven useful for many atmospheric applications. Some of these methods include the Kalman filter (KF) for linear dynamics and the extended Kalman filter (EKF) for nonlinear dynamics, and to improve the efficiency of the traditional EKF, another method called the ensemble Kalman filter (EnKF) was developed. The 4DVAR and 3DVAR algorithms are the two main variational DA methods that tend to combine the observations and numerical models in an optimal way to give the best possible estimate of the initial model state by minimizing the cost function that measures the weighted sum of square errors. When the time-dependent dimension is taken out, 4DVAR becomes 3DVAR making the former a better scheme for dynamical systems than the latter. It was also revealed that although variational schemes are rarely used, they are found to be more successful than other algorithms (Ide et al.,1997). Particularly, the 4DVAR data assimilation method that seeks to minimize the cost function defined as the weighted sum of squared deviation between the model and the observations became operational in 1997 Eknes and Evensen (1997). This method corrects the model in a fourdimensional domain, that is, the three spatial dimensions and the time dimension using the model dynamics as a constraint, either strong or weak. The two approaches to this method are the weak and strong constraints 4DVAR and these differ in such a way that the strong constraints 4DVAR assume a perfect model implying the error is confined in the initial model state while the weak constraints assume errors in the observations, model external forcing, initial and boundary conditions relaxing the perfect model assumption as discussed by Gustafsson et al. (2012). According to Zupanski (1997), the weak constraints 4DVAR has always yielded accurate results because it includes all the information that may be lost by the strong constraints method and its also capable of altering the dynamics of a given model hence fitting the observations. The effectiveness of weak constraints 4DVAR algorithm in correcting and estimating the model error in 4DVAR for long-range assimilation windows was the motivation for its use in this study as it assimilates observations through the minimization of a least-squares cost function which is a function of the error between the model and the measurements, seeking for an optimal solution. Section 2.2. Overview of weak constraint 4D-Var DA technique Page 5 To solve the dynamic problem, coupled Euler Lagrange equations are obtained from the minimized problem for which a convergent iteration can then be defined. The improved estimates of the poorly known parameters can be calculated by solving the problem for each of the linear iterates using the representer which expresses the problem as the sum of a first guess and a finite linear combination of representer functions. It formulates the problem in such a way that the solution procedure involves only covariance multiplications, and not their inverses, with variables defined in the state and the observations spaces. The representer method is adopted for solving weak constraints 4DVAR problems because of its ability to decouple the Euler Lagrange system hence giving the minimum solution as portrayed by Bennett (1992). 3. Weak Constraint 4DVAR This chapter discusses the implementation of weak constraint 4DVAR with representers to assimilate wind fields, emmision, and PM2.5concentrations data directly into the transport model (1.2.1) as it allows to alter the dynamics of the model for better results. Using this method involved derivation of Euler Lagrange equations when optimizing the state variable q(x, t)as well as the input that is; Q(x, t),u(x, t),K, and cand then using the representer method to decouple the E-L systems making the scheme computationally feasible. 3.1 The model In this study, the model used for estimating P M2.5concentrations q(x, t) is an Eulerian one dimensional transport model given as below; (qt(x, t) + ux(x, t)q(x, t) + u(x, t)qx(x, t)−Kqxx(x, t) + cq(x, t) = Q(x, t) q(x, 0) = q0(x)x∈Ω, t ∈[0, T](3.1.1) where q=q(x, t)is the concentration of PM2.5in (µgm−3), q0(x,t) is the initial concentration of PM2.5in (µ gm−3), Q=Q(x, t)is the source term which denotes the rate of mass influx per unit volume per unit time into the smoke plume (kgm−3s−1), u=u(x, t)represents dispersion by ambient wind in (ms−1), Kis a one dimensional diffusivity constant representing turbulent diffusion in (m2s−1) and cis the decay coefficient of the fall out due to rain or temperature in (Bq)and since this is a one dimensional model. If inputs u(x, t),K,cand Q(x, t)are represented accurately, equation (3.1.1) contains all of the factors that contribute to smoke behavior. This includes, the emission source, plume rise, transport and dispersion by ambient wind, as well as chemical reactions in the smoke. 3.2 Data The data assumes a finite number of observations collected within a temporal domain 0≤t≤T and a spatial domain 0≤x≤L. Following Furtado et al. (2010), these observations are imperfect point measurements of the independent variable β(x, t)collected at M points in space at time (xm, tm)and these are given by; dm=Hββ(xm, tm) + m1≤m≤M where Hβis the observation operator that transforms the unknown β(x, t)into observation equivalents and mis the measurement error. 6 Section 3.4. Derivation of Euler Lagrange systems separately Page 13 where λ(x, t) = Wf∂q ∂t +∂u ∂x q+u∂q ∂x −K∂2q ∂x2−c∂q ∂x −ˆ Q(x, t). The filtering property of the Dirac delta was then employed to equation (3.4.15) to eliminate δ(xj, tj); J X j=1 (dQ)j−HQδQ(xj, tj)ˆ Q(xj, tj)HT Q=ZT 0ZΩ J X j=1 (dQ)j−HQˆ Q(xj, tj) HT Qδ(x−xj)δ(t−tj)δQdxdt (3.4.16) Having expressed δJin terms of δQ, and as argued that δJ=O(δQ)2, for ˆ Qto be the local extremum of J, the coefficients of δQ(x,t) must vanish giving the E-L equations below;        ˆ Q(x, t) = Q0(x, t) + CEλ(x, t) + CEWQPJ j=1 (dQ)j−HQˆ Q(xj, tj)HT Qδ(x−xj)δ(t−tj) λ(x, t) = Wf∂q ∂t +∂u ∂x q+u∂q ∂x −K∂2q ∂x2−c∂q ∂x −ˆ Q(x, t)(3.4.17) Solving the E-L system The optimal estimate ˆ Q(x, t)could then be obtained from the E-L equations by substituting forλ(x, t)in the equation for ˆ Qabove as shown below: ˆ Q(x, t) =                  Q0(x, t) + CEWf∂q ∂t +∂u ∂x q+u∂q ∂x −K∂2q ∂x2−c∂q ∂x −ˆ Q(x, t);for x 6=xj Q0(x, t) + CEWf∂q ∂t +∂u ∂x q+u∂q ∂x −K∂2q ∂x2−c∂q ∂x −ˆ Q(x, t) +CEWQ((dQ)j−HQˆ Q(xj, tj); for x =xj (3.4.18) Simplifying equation (3.4.18) gives; ˆ Q(x, t) =            A−1Q0(x, t) + CEWf∂q ∂t +∂u ∂x q+u∂q ∂x −K∂2q ∂x2−c∂q ∂x ;for x 6=xj B−1Q0(x, t) + CEWf∂q ∂t +∂u ∂x q+u∂q ∂x −K∂2q ∂x2−c∂q ∂x +CEWQdQ;for x =xj (3.4.19) where A=I+CEWf,B=I+CEWf+CEWQHQand Iis an identity matrix 3.4.3 Optimizing for the diffusivity constant. In this sub-section, the E-L system was derived by minimizing the cost function (3.3.4) with respect to K. Similarly the calculus of variations technique was employed as earlier to minimize the functional in a least square sense hose the local extremum was assumed to be ˆ Khence giving the first variation δJ=J[ˆ K+δK]− J [ˆ K] Section 3.4. Derivation of Euler Lagrange systems separately Page 14 as below: δJ=−2WfZT 0ZΩ∂q ∂t +∂u ∂xq+u∂q ∂x −ˆ K∂2q ∂x2−c∂q ∂x −Q∗∂2q ∂x2δKdxdt +2WΓ(ˆ K−K0)δK (3.4.20) let λ(x, t) = WfRT 0RΩ∂q ∂t +∂u ∂x q+u∂q ∂x −ˆ K∂2q ∂x2−c∂q ∂x −Q∗∂2q ∂x2dxdt. As argued by Bennett (2005) that δJ=O(δK)2, for ˆ Kto be a local extremum of Jthen coefficients of δK in equation (3.4.20) must go to zero giving the E-L system below;    ˆ K=K0+CΓλ(x, t) λ(x, t) = WfRT 0RΩ∂q ∂t +∂u ∂x q+u∂q ∂x −ˆ K∂2q ∂x2−c∂q ∂x −Q∗∂2q ∂x2dxdt (3.4.21) Solving the E-L system The optimal estimate ˆ Kwas obtained by substituting for λin the equation for ˆ Kas shown below ˆ K=K0+CΓWfZT 0ZΩ∂q ∂t +∂u ∂xq+u∂q ∂x −ˆ K∂2q ∂x2−c∂q ∂x −Q×∂2q ∂x2dxdt (3.4.22) Making ˆ Kthe subject in equation (3.4.22) above yields: ˆ K= K0+CΓWfRT 0RΩ∂q ∂t +∂u ∂x q+u∂q ∂x −c∂q ∂x −Q×∂2q ∂x2dxdt 1 + CΓWfRT 0RΩ ∂2q ∂x2×∂2q ∂x2dxdt (3.4.23) 3.4.4 Optimizing for the decay coefficient. For this case, the cost function was minimized in the least square sense with respect to cusing the calculus of variations technique seen earlier to obtain the maximum likelihood estimate for the decay coefficient. The local extremum for the functional (3.3.4) was assumed to be ˆcand hence the first variation δJ=J[ˆc+δc]− J [ˆc]was obtained as; δJ=−2WfZT 0ZΩ∂q ∂t +∂u ∂xq+u∂q ∂x −K∂2q ∂x2−ˆc∂q ∂x −Q×∂q ∂xδcdxdt + 2Wγ(ˆc−c0)δc (3.4.24) let λ(x, t) = WfRT 0RΩ∂q ∂t +∂u ∂x q+u∂q ∂x −K∂2q ∂x2−ˆc∂q ∂x −Q×∂q ∂x dxdt. Basing on the argument that δJ=O(δc)2, for ˆcto be a local extremum of Jthen coefficients of δc in equation (3.4.24) must vanish giving the E-L system of equations as;    ˆc=c0+Cγλ(x, t) λ(x, t)WfZT 0ZΩ∂q ∂t +∂u ∂xq+u∂q ∂x −K∂2q ∂x2−ˆc∂q ∂x −Q×∂q ∂xδcdxdt (3.4.25) Section 3.4. Derivation of Euler Lagrange systems separately Page 15 Solving the E-L system Substituting for λ(x, t)in the equation for ˆcand then making ˆcthe subject gives: ˆc= c0+CγWfRT 0RΩ∂q ∂t +∂u ∂x q+u∂q ∂x −K∂2q ∂x2− −Q×∂q ∂x dxdt 1 + CγWfRT 0RΩ ∂q ∂x ×∂q ∂x dxdt (3.4.26) 3.4.5 Optimizing for the dispersion by ambient wind. For this case the E-L system of equations was derived when optimizing for the dispersion by ambient wind and this was done by minimizing the cost function (3.3.4) in the least square sense with respect to u(x, t)to obtain the maximum likelihood estimate. The calculus of variations method was used following similar steps as those for obtaining equation (3.4.3) above with ˆu= ˆu(x, t)as the local extremum to the functional (3.3.4) above giving δJ=J[ˆu+δu]− J [ˆu]as below; δJ= 2 ZT 0ZΩ λ(x, t)∂δu ∂x q(x, t)dxdt + 2 ZT 0ZΩ λ(x, t)qxδudxdt + 2WEuZT 0ZΩˆu(x, t)−u0(x, t)δudxdt −2Wu N X n=1 (du)n−Huˆu(xn, tn)HT uδu (3.4.27) where λ(x, t)≡Wf∂q ∂t +∂ˆu ∂x q+ ˆu∂q ∂x −K∂2q ∂x2−c∂q ∂x −Q(x, t). Unlike equation (3.4.15) whose terms were only in δQ, for (3.4.27) the terms are both in δu and ∂δu ∂x hence the need for integrating the equation by parts and also using the filtering property of the Dirac delta functions to express δJin terms of δu only as shown below: δJ= 2 ZT 0 λ(L, t)q(L, t)δu(L, t)dt −2ZT 0 λ(0, t)q(0, t)δu(0, t)dt −2ZT 0ZΩ ∂λ ∂xq(x, t)δudxdt + 2WEuZT 0ZΩˆu−u0(x, t)δudxdt −2WuZT 0ZΩ N X n=1 (du)n−Huˆu(xn, tn)HT uδ(x−xn)δ(t−tn)δudxdt (3.4.28) The boundary terms are equal and opposite, due to the periodic boundary condition. Having expressed δJentirely in terms of δu and as argued earlier that δJ=O(δu)2, this implies that ˆu is the extremum of Jprovided the coefficients of δu(x,t), δu(0,t) and δu(L,t) vanish hence giving the EulerLagrange equations below;              −∂λ ∂x q(x, t) = −WEu[ˆu(x, t)−u0(x, t)] + WuPN n=1 (du)n−Huˆu(xn, tn)HT uδ(x−xn)δ(t−tn) λ(L, t) = λ(0, t) ∂ˆu ∂x q+ ˆu∂q ∂x =Cfλ(x, t)−∂q ∂t +K∂2q ∂x2+c∂q ∂x −Q(x, t) q(L, t) = q(0, t) (3.4.29) Section 3.5. Derivation of the Euler Lagrange system simultaneously Page 16 Solving the E-L system (B)   −∂λ ∂x q(x, t) = −WEu[ˆu(x, t)−u0(x, t)] + WuPN n=1 (du)n−Huˆu(xn, tn)HT uδ(x−xn)δ(t−tn) λ(L, t) = λ(0, t) (3.4.30) (F)(∂ˆu ∂x q+ ˆu∂q ∂x =Cfλ(x, t)−∂q ∂t +K∂2q ∂x2+c∂q ∂x −Q(x, t) q(L, t) = q(0, t)(3.4.31) This is solved using numerical finite difference schemes let λk i≈λ(xi, tk),qk i≈q(xj, tk), and ˆuk i≈ˆu(xj, tk)for 1≤i≤Iand 1≤k≤K the solution can be obtained from: λk i+1 −λk i ∆xqk i=(−WEu[ˆuk i−u0k i]for xi6=xn −WEu[ˆuk i−u0k i] + Wudu−Huˆuk ifor xi6=xn (3.4.32) ˆuk i+1 −ˆuk i ∆xqk i=−+qk i−qk i−1 ∆xuk iCfλk i−fk i(3.4.33) 3.5 Derivation of the Euler Lagrange system simultaneously This section shows the derivation of the E-L system when optimizing for the concentration of PM2.5, emission source, ambient dispersion of wind, the decay constant due to rainfall or temperature and the diffusivity constant simultaneously. This was done by minimizing the cost function (3.3.4) in a least square sense with respect to q(x, t),Q(x, t),u(x, t),Kand cusing the calculus of variations technique as earlier used when optimizing for q(x, t). For this, the local extrema for the functional Jwere assumed to be ˆq,ˆ Q,ˆu,ˆ Kand ˆcrespectively. The first variation; δJ =J[ˆq+δq, ˆ Q+δQ, ˆu+δu, ˆ K+δK, ˆc+δc]−J[ˆq, ˆ Q, ˆu, ˆ K, ˆc]was then obtained as Section 3.5. Derivation of the Euler Lagrange system simultaneously Page 17 below; δJ= 2WfZT 0ZΩ∂ˆq ∂t +∂ˆu ∂x ˆq+ ˆu∂ˆq ∂x −ˆ K∂2ˆq ∂x2−ˆc∂ˆq ∂x −ˆ Q× ∂δq ∂t +∂ˆu ∂xδq +∂δu ∂x ˆq+δu ∂ˆq ∂x + ˆu∂δq ∂x −δK ∂2ˆq ∂x2−ˆ K∂2δq ∂x2−δcˆqx−ˆc∂δq ∂x −δQdxdt + 2WEZT 0ZΩˆ Q−Q0(x, t)δQdxdt + 2WEuZT 0ZΩˆu−u0(x, t)δudxdt + 2WiZΩˆq(x, 0) −q0(x, t)δqdx −2Wq M X m=1 (dq)m−Hqˆq(xm, tm)HT qδq− 2WQ J X j=1 (dQ)j−HQˆ Q(xj, tj)HT QδQ −2Wu N X n=1 (du)n−Huˆu(xn, tn)HT uδu + 2WΓ(ˆ K−K0)δK + 2Wγ(ˆc−c0)δc (3.5.1) let λ=λ(x, t)≡Wfˆqt(x, t) + ˆuxˆq+ ˆuˆqx−ˆ K∇2ˆq(x, t)−ˆc∇ · ˆq(x, t)−ˆ Q. Since equation (3.5.1) is in terms of δq,δu,δQ,δK,δc,∂δq ∂t ,∂δq ∂x , and ∂δu ∂x , integration by parts as well as the filtering property of the Dirac delta functions were employed to make the equation be in terms of δq,δu,δQ,δK, and δc only as shown below; Section 3.5. Derivation of the Euler Lagrange system simultaneously Page 18 δJ= 2 ZΩλδqT 0 dx −2ZT 0ZΩ ∂λ ∂t δqdxdt + 2 ZT 0ZΩ λ∂ˆu ∂xδqdxdt + 2 ZT 0λˆqδuΩ dt −2ZT 0ZΩλ∂ˆq ∂x +∂λ ∂x ˆqδudxdt + 2 ZT 0λˆuδqΩ dt −2ZT 0ZΩλ∂ˆu ∂x +∂λ ∂x ˆuδqdxdt + 2 ZT 0ZΩ λ∂ˆq ∂xδudxdt −2ˆ KZT 0λ∂δq ∂x Ω dt + 2 ˆ KZT 0∂λ ∂xδqΩ dt −2ˆ KZT 0ZΩ ∂2λ ∂x2δqdxdt −2ZT 0ZΩ λ∂2ˆq ∂x2δKdxdt −2ˆcZT 0λδqΩ dt + 2ˆcZT 0ZΩ ∂λ ∂xδqdxdt −2ZT 0ZΩ λ∂ˆq ∂xδcdxdt −2ZT 0ZΩ λδQdxdt + 2WEZT 0ZΩˆ Q−Q0(x, t)δQdxdt + 2WEuZT 0ZΩˆu−u0(x, t)δudxdt + 2WiZΩˆq(x, 0) −q0(x, t)δqdx −2WqZT 0ZΩ M X m=1 (dq)m−Hqˆq(xm, tm)HT qδ(x−xm)δ(t−tm)δqdxdtx −2WQZT 0ZΩ J X j=1 (dQ)j−HQˆ Q(xj, tj)HT Qδ(x−xj)δ(t−tj)δQdxdt −2WuZT 0ZΩ N X n=1 (du)n−Huˆu(xn, tn)HT uδ(x−xn)δ(t−tn)δudxdt +WΓ(ˆ K−K0)δK +Wγ(ˆc−c0)δc (3.5.2) Taking into account that the boundary terms are equal and opposite, due to the periodic boundary condition and the argument made by Bennett (2005) that δJ=O(δq)2+O(δQ)2+O(δu)2+ O(δK)2+O(δc)2, for ˆq,ˆ Q,ˆu,ˆ Kand ˆcto be local extrema for J, coefficients of δq(x, t),δq(x, 0), δq(x, T),δq(0, t),δq(L, t),δu(x, t),δu(0, t),δu(x, L),δQ(x, t),δK(x, t), and δc(x, t)must vanish giving the system of Euler Lagrange equations below: Section 3.5. Derivation of the Euler Lagrange system simultaneously Page 19                                                      −∂λ ∂t −ˆ K∂2λ ∂x2−(ˆu(x, t)−ˆc)∂λ ∂x =WqPM m=1 (dq)m−Hqˆq(xm, tm)HT qδ(x−xm)δ(t−tm) ∂ˆq ∂t +∂ˆu ∂t ˆq+ ˆu∂ˆq ∂t −ˆ K∂2ˆq ∂x2−ˆc∂ˆq ∂t =ˆ Q+Cfλ(x, t) ˆu(x, t) = u0(x, t) + CEu ∂λ ∂x ˆq(x, t) + CEuWuPN n=1 (du)n−Huˆu(xn, tn)HT uδ(x−xn)δ(t−tn) ˆ Q(x, t) = Q0(x, t) + CEλ(x, t) + CEWQPJ j=1 (dQ)j−HQˆ Q(xj, tj)HT Qδ(x−xj)δ(t−tj) ˆ K=K0+CΓRT 0RΩλ(x, t)∂2ˆq ∂x2dxdt ˆc=c0+CγRT 0RΩλ(x, t)∂ˆq ∂x dxdt ˆq(x, 0) = q0(x, t) + Ciλ(x, 0) λ(x, T) = 0 λ(L, t) = λ(0, t) ˆq(L, t) = ˆq(0, t) (3.5.3) Solving the E-L system Re ordering the EulerLagrange of equation (3.5.3) gives the backward problem (3.5.4) , forward problem (3.5.5) and the optimal estimates of the parameters as shown below: (B)         −∂λ ∂t −ˆ K∂2λ ∂x2−(ˆu(x, t)−ˆc)∂λ ∂x =WqPM m=1 (dq)m−Hqˆq(xm, tm)HT qδ(x−xm)δ(t−tm) λ(x, T) = 0 λ(L, t) = λ(0, t) (3.5.4) (F)     ˆqt(x, t) + ˆux(x, t)ˆq+ ˆuˆqx−ˆ K∇2ˆq(x, t)−ˆc∇ · ˆq(x, t) = ˆ Q+Cfλ(x, t) ˆq(x, 0) = q0(x, t) + Ciλ(x, 0) ˆq(0, t) = ˆq(L, t) (3.5.5) with parameter estimates;                  ˆu(x, t) = u0(x, t) + CEu ∂λ ∂x ˆq(x, t) + CEuWuPN n=1 (du)n−Huˆu(xn, tn)HT uδ(x−xn)δ(t−tn) ˆ Q(x, t) = Q0(x, t) + CEλ(x, t) + CEWQPJ j=1 (dQ)j−HQˆ Q(xj, tj)HT Qδ(x−xj)δ(t−tj) ˆ K=K0+CΓRT 0RΩλ(x, t)∂2ˆq ∂x2dxdt ˆc=c0+CγRT 0RΩλ(x, t)∂ˆq ∂x dxdt (3.5.6) 4. Pseudo Code for the E-L Systems 4.1 Sample pseudo code for ˆq Basing equations (3.4.8) to (3.4.14) above, the pseudo code for ˆqcan be represented as below: 1. Start 2. Input covariance matrix Cqand matrix for the observation operator Hq 3. Compute qFby numerical integration of equation (3.4.9) 4. Compute innovation haccording to h=dq−HqqF 5. Compute the adjoint representer αm(x, t)according to equation(3.4.10) 6. Compute the representer rm(x, t)for 1≤m≤Maccording to equation(3.4.11) 7. compute rm(xl, tl)for m= 1,2,· · · , M and l= 1,2,· · · , m that is matrix RM×M 8. Compute P= (R+CqI) 9. Compute the analysis given by the following equation: ˆq(x, t) = qF(x, t) + hTP−1r(x, t) 10. Stop 4.2 Sample pseudo code for ˆ Q ˆ Q(x, t) =            A−1Q0(x, t) + CEWf∂q ∂t +∂u ∂x q+u∂q ∂x −K∂2q ∂x2−c∂q ∂x ;for x 6=xj B−1Q0(x, t) + CEWf∂q ∂t +∂u ∂x q+u∂q ∂x −K∂2q ∂x2−c∂q ∂x +CEWQdQ;for x =xj (4.2.1) where A=I+CEWf,B=I+CEWf+CEWQHQand Iis an identity matrix The pseudo code for (4.2.1) is as below: 1. Start 2. Define Function Q0(x,t) 3. Define Function Q(x,t) 20 Section 4.2. Sample pseudo code for ˆ QPage 21 4. Input dQ, xm,CE,Wf, and WQ 5. Input lower limit (0), upper limit (L) and number of sub intervals (N) for x 6. Input lower limit (0), upper limit (T) and number of sub intervals (M) for t 7. Calculate A and B 8. Create arrays for t and x using their lower limits, upper limits and number of sub intervals respectively 9. Create an empty array for ˆ Q(x, t)with dimension (M+1)×(N+ 1) 10. Start a for loop as shown below: for t= 0 toTdo for x= 0 to Ldo if x6=xjthen ˆ Q(x, t) = A−1Q0(x, t) + CEWfQ(x, t) elseˆ Q(x, t) = B−1Q0(x, t) + CEWfQ(x, t) + CEWQdQ end if end for end for 11. Print ˆ Q(x, t)as the result 12. Stop The code above was implemented to simulate the emission source as shown in figure 4.1 below. The covariance matrices in (4.2.1) were taken to be equal to one, the initial source term was assumed to be an exponential term since smoke grows/ decays exponentially that is: Q0(x, t) = cexp(x−xj) ∂q ∂t +∂u ∂xq+u∂q ∂x −K∂2q ∂x2−c∂q ∂x =c1exp(x−xj) and data points dQwere also assumed, c= 1 and c1= 0.96. Section 4.3. Sample pseudo code for ˆ KPage 22 Figure 4.1: Plot of exact and estimate emission sources at different points in space. 4.3 Sample pseudo code for ˆ K ˆ K= K0+CΓWfRT 0RΩ∂q ∂t +∂u ∂x q+u∂q ∂x −c∂q ∂x −Q×∂2q ∂x2dxdt 1 + CΓWfRT 0RΩ ∂2q ∂x2×∂2q ∂x2dxdt (4.3.1) Basing equation (4.3.1) above, the pseudo code for ˆ Kcan be represented as below: 1. Start 2. Define Function f1(x, t) = qt(x, t) + ∇ · u(x, t)q(x, t)−c∇ · q(x, t)−Q(x, t)*∇2q(x, t) 3. Define Function f2(x, t) = ∇2q(x, t)∗ ∇2q(x, t) 4. Input K0,CΓ, and Wf 5. Input lower limit (0), upper limit (L) and number of sub intervals (N) for x 6. Input lower limit (0), upper limit (T) and number of sub intervals (M) for t 7. Create arrays for t and x using their lower limits, upper limits and number of sub intervals respectively 8. To calculate for the double integral in the numerator in equation (4.3.1)the Trapezoidal rule can be used to give an output Result 1 9. To calculate for the double integral in the denominator in equation (4.3.1)the Trapezoidal rule can be used to give an output Result 2