Most of our knowledge of flow in porous media is obtained at the pore and the macro scale. For reservoir scale modelling it is not practical to model the flow at these fine scales. Considering the usual objectives (e.g. large scale flow pattern and production forecast) it is undesirable to have to gather the tremendous amount of fine scale data that is required to model an entire reservoir. Moreover, present computational resources are simply not able to handle flow simulations of this size. Hence we resort to models that describe the essential physical behaviour in an averaged sense at the mega scale without modelling all finer scale details. Unfortunately, fine scale details are often correlated over large distances and appear to be part of the essential behaviour at a larger scale. As a consequence the reservoir scale flow is not just a function of some average properties of the finer scale structures. For these cases the scales in the problem are non-separable. The major challenge in upscaling is to account for this non-separability of scales. Going from single phase single component, via single phase multi component to multi phase fluid flow, upscaling becomes progressively more complex. In this Thesis we treat several aspects of upscaling of single and multi phase flow in porous media approached from a single phase context. Single phase upscaling The simplest flow type in porous media is single phase single component flow. For this type of flow the general assumption is that the upscaled megascopic pressure flux relation is the same as at the macro scale, namely Darcyâs law, albeit with effective coefficients. When the spatial length scales are non-separable the upscaled quantity depends on choices we make at the fine scale. As a consequence, single phase upscaling is mainly concerned with the proper estimation of the fine scale flux boundaries. In this thesis a method is developed that allows us to investigate the pressure and flux response of each individual flux boundary independent of the other boundaries. Fourier decomposition enables us further to investigate the effect of each individual frequency a flux boundary consists of. The major observation is that the significance of the flux response in the far field generally decreases with increasing frequency of the flux boundary that is applied. This does not generally hold for the pressure response; In case of long correlated streaks of especially low permeability the far field pressure response may be significantly larger than one would a priori expect. Upscaling of tracer dispersion At the next level of complexity we arrive at the upscaling of tracer dispersion. Traditionally dispersion has been modelled using Fickian models. However the Fickian approach does not satisfactorily explain all of the experimental observations. For instance the Fickian model cannot describe effective dispersion coefficients that vary with the Péclet number or with the length scale of observation. Moreover the Fickian description cannot account for the partial reversibility of dispersion when the flow direction is reversed. We focus on (Taylor) dispersion of a tracer released in a unidirectional velocity field belonging to a two dimensional (2D) porous medium. The tracer concentration is described by the two dimensional unidirectional classical convection dispersion equation (2D uCCDE), in which we initially restrict the small scale dispersive mechanisms to isotropic molecular diffusion. Our goal is to describe the (non-Fickian) evolution of the height averaged tracer concentration in time. The analysis starts with the relaxation concept; transverse diffusion causes the particle cloud to describe a transition from a correlated convective behaviour for short times towards uncorrelated Fickian behaviour for asymptotic long times. From a particle viewpoint this process describes the convergence of the velocity auto correlation function (VACF) of tracer particles to zero in time, as a result of the transverse displacement by diffusion. The correlation of the velocity of a particle in time results from the spatial correlation of the velocity field in the transverse direction and the limited (averaged) transverse displacement of particles by diffusion over short times. This relaxation process is characterised by a height averaged concentration profile that initially reflects the velocity profile and that converges in time towards a Gaussian profile typical for Fickian displacements. The corresponding spatial variance is initially proportional to t2 and becomes asymptotically with time proportional to t. The time that characterises this relaxation process is called the relaxation time. To analyse the model behaviour we transform the 2D uCCDE into an equivalent spectral representation using Fourier transformation. Each non-zero mode contributes to the mass balance of the height averaged concentration (zero-th mode) via a dispersive modal flux term. The evolution of each non-zero mode is described by a relaxation equation. It is characterised by a modal relaxation time that originates from the transverse molecular diffusion term in the 2D uCCDE. This modal relaxation time is proportional to the total height of the velocity field squared and decreases inversely proportional to the molecular diffusion and to the modal number squared. From the spectral representation we are able to derive the exact behaviour of the spatial moments belonging to the height averaged concentration. Moreover it forms the base of the upscaling approach we developed. To obtain the essential behaviour of the height averaged concentration, it is impractical to solve all modal equations. Hence we need to reduce the set of equations to a manageable size. A common approach is to describe the evolution of the dispersive Taylor flux as a single valued quantity. An important consequence of this approach is that the multi-scale character of the full model is lost. Classically this Taylor flux has been modelled in a Fickian way. However, such an approach ignores the relaxation process of the tracer and only covers the tracer behaviour for asymptotic long times. Here we follow the upscaling approach of Camacho. By summing all evolution equations of the modal fluxes, an approximate evolution equation of the Taylor flux is obtained. The result is a linear parabolic relaxation equation that is characterised by an approximate effective relaxation time and a closure term that accounts for all higher order modal interactions. Combining this equation with the mass balance for the averaged concentration results in the so-called c-J model. In case the spreading by diffusion in the longitudinal direction is negligible compared to the spreading induced by the velocity field, the c-J model simplifies to a generalised hyperbolic Telegraph equation. The c-J model describes the variance correctly in both the short and the long time limit. As it incorporates the relaxation process in a single scale sense, it is also able to show qualitatively the proper development of the variance for intermediate times. Quantitatively, we show that the c-J model improves significantly when scale separation is applied; The smaller Fourier scales that are relaxed with respect to the time scale of observation are described by a Fickian term and only the relaxation of the large(r) scales is modelled explicitly. In this way part of the multi scale character of the full model is again retrieved. Special attention is paid to the closure term in the c-J model that accounts for the modal interactions of the higher order concentration modes. As mentioned the multi scale character of the 2D uCCDE is lost by describing the dispersive flux as a single valued quantity. A consequence of this is that the c-J model is unable to show a sign change in the 3-rd moment that may occur in the full 2D model. To remedy this we consider alternative descriptions of the closure term in the upscaled model. We restrict our focus to the proper description of the third spatial moment and require that the present description of the mean and variance remains unaltered. Moment analysis demonstrates that the present form of the term that accounts for the modal interactions in the c-J modal is the only allowable linear form that does not affect the present description of the mean and variance. Alternatively, the coefficient in front of the closure term may be replaced by an empirical time varying function. The result describes qualitatively the proper behaviour of the 3-rd moment and is exact in the convective and dispersive limit. Overall it produces smaller errors in the 3-rd moment than the original c-J model while the lower moments remain unchanged. Until now we discussed the behaviour for particle distributions that are initially uniform over the height. Distributions that are initially non-uniform converge asymptotically in time towards a uniform distribution. This is shown by the convergence of the mean particle velocity towards the mean fluid velocity in time. Another consequence of non-uniform distributions is that the modal interactions of the higher order concentrations modes also affect the variance. Moreover, each moment is characterised by its "own" effective relaxation time. In general it is not longer adequate to describe the behaviour of the modal fluxes as a single valued quantity. We develop a truncation model in which, up to a separation scale, the evolution of each mode is modelled individually while the smaller scales are combined in an effective upper mode. The separation scale is determined such that the mean and variance are correctly described. The detail required in the truncation model increases as the initial distribution deviates more from uniformity and increases for decreasing observation times. However, even for quite small observation times, the truncation model is able to describe the proper behaviour with only little detail. We compare the 1D CCDE and the Telegraph model with the truncation method by means of the first two spatial moments. The CCDE ignores the relaxation process and only describes the behaviour of the mean and variance properly for asymptotic long times. For non-uniform initial particle distributions, the Telegraph equation describes both a relaxation of the mean and the variance with a single effective relaxation time. It is exact in the convective limit and describes the proper growth of the mean and variance in the long time limit. However, the shortcoming of the Telegraph model is that it cannot distinguish two distributions that have the same mean particle velocity and, as a consequence, it may produce large errors for intermediate times. Its applicability is therefore, similarly to the 1D CCDE, essentially restricted to uniform initial particle distributions. A second concept that is incorporated in the full physical model is reversibility of dispersion. When the flow direction is reversed a 'new' relaxation process is initiated similar to forward relaxation. However, the velocity of a particle after flow reversal is correlated in time to the velocity of that particle before reversal. This is again caused by the limited transverse particle displacement by diffusion over short times and the correlation of the velocity field in the transverse direction. As a consequence, we experience a demixing of the particle cloud as long as this velocity correlation dominates the particle behaviour. This is quantified by a temporary decrease of the variance. The major advantage of the c-J model over the 1D Fickian CCDE is demonstrated in these cases of partial reversibility of dispersion. While the Fickian model simply cannot handle a temporal decrease of the variance, the reversal behaviour is directly incorporated in c-J model via the relaxation time. The c-J model describes the (reversed flow) variance exactly in the convective and dispersive limit and gives a qualitatively correct picture for intermediate times. Applying scale separation further improves the quantitative description for intermediate times. So far the small scale mixing mechanism was restricted to isotropic diffusion. For flow in porous media this is too restrictive and we generalise it to a transverse varying dispersion. For unidirectional flow the microscopic dispersion in the transverse direction is commonly defined in an averaged sense at the macro scale as . This has as consequence that the relaxation mechanism turns from a time dependent mechanism for small Péclet number towards a space dependent relaxation mechanism for large Péclet number. For the spectral representation the extended dispersion has as major consequence that the relaxation times of the higher order modes are coupled. However, we have shown that by transformation of the modal concentrations we can obtain a set of modified modal equations that is of identical shape as the original set in which the modified relaxation times are again uncoupled. The transformation depends on the Péclet number. The advantage of the shape similarity is that, except for the modified coefficients, the mean and variance for uniform particle distributions are not altered and only one additional term shows up in the variance for non-uniform particle distributions. An upscale approach is followed similar to the diffusive case. The result is again a parabolic relaxation equation for the Taylor flux that is similar to the case of isotropic molecular diffusion. However, the effective relaxation time and modal interaction term are now functionally related to the Péclet number. In the diffusive limit the upscaled equation describes a relaxation in time identical to the isotropic case. For large Péclet numbers it describes a relaxation in space, similar to the full 2D model. We also investigate tracer dispersion in arbitrary velocity fields belonging to two dimensional porous media bounded in the y-direction. In these fields a second convergence mechanism may be present of which the variance of the height averaged concentration is initially proportional to t2 and which is linear in time for asymptotic long times. This mechanism, called convective convergence, takes place if the velocity field is second order stationary and if the velocity auto correlation function converges to zero sufficiently fast in time. As this mechanism is purely convective its dispersion is in principle fully reversible. In other words the spreading caused by the forward movement is completely cancelled when the flow is reversed and the flow time in both directions is equal. We investigate again if it is possible to describe the dispersive flux as a single valued quantity. Starting from the spectral equivalent of the two-dimensional convection diffusion equation, we follow an upscaling approach similar to the one discussed above. However, the spatial dependency of the velocity complicates matters considerably. Even for particle distributions that are initially uniform over the height, the modal interactions affect the lowest order derivatives in the equations of the dispersive flux and as a consequence the complete relaxation mechanism. The resulting model is of identical form as the c-J model but with coefficients that vary in space and are all functions of the velocity. If the velocity field is unidirectional the original c-J model is retrieved. In the absence of diffusion the equation is purely convective and is from a one-dimensional viewpoint able to demonstrate convective convergence and its reversibility. However, the present model has some severe (physical) shortcomings. Convective convergence in the sense of a single valued dispersive flux, requires that the relaxation time alters sign when the flow is reversed. A negative relaxation time is, similar to a negative dispersion coefficient hard to defend. Besides the relation between the upscaled coefficients and the original two dimensional field is (yet) unclear. Nevertheless, the present extension to arbitrary velocity field appears as a promising ansatz. The megascopic equations that we derive all appear to be of the same form, the c-J model, with the exception of the truncation model. The c-J model is also similar to recently developed macroscopic non-Fickian dispersion models. However the c-J model is a 1D model whereas the others are 3D. This similarity encouraged us to investigate if experimental results for macroscopic homogeneous media can be explained with this model. When a relaxation mechanism is present that works similarly to the combination of dispersion and diffusion, we show indeed that the c-J model may explain the variation of the effective dispersion coefficient as function of the Péclet number. It may also explain some experimental results concerning the reversibility of dispersion. However at present a clear mathematical link between the 1D c-J model and macroscopic homogeneous media is not yet established. Finally we establish a link between the relaxation process in the 2D uCCDE and Markovian Lagrangian theory. A Markovian model is presented in which we explicitly modelled the loss of velocity correlation of particles in time. A discrete velocity set represents the unidirectional velocity field of the full two-dimensional model. A large number of particles are redistributed over this set in time. If the magnitude of the small scale dispersion mechanisms in the full 2D model decreases the redistribution is such that the probability for a particle to keep the same velocity (auto correlation) increases. The model is constructed in such a way that any arbitrary initial particle distribution converges towards a uniform particle distribution in time. When the velocity auto correlation is uniform, the model shows relaxation behaviour that is essentially the same as the Telegraph equation. Moreover, in case there are only two velocities we proof that both models are identical upto second order accuracy in time. For a non-uniform velocity auto correlation, multi scale behaviour is demonstrated that is similar to the 2D uCCDE. By exploring a moment analysis from a Lagrangian viewpoint we derive analytical expressions for the mean, variance and third moment. More importantly, the Lagrangian viewpoint also enables us to explain the behaviour of the moments of the 2D uCCDE qualitatively. Especially the correlation structure of the velocity in time is visible from this point of view. Two phase flow upscaling This dissertation ends with upscaling of the more complex two phase flow type from a single phase context. In two phase flow in heterogeneous porous media, the important âconvectiveâ part may be divided in a single phase term and a two phase flow term. The first term captures the large scale permeability variations and is expressed in a Piston-Like (single phase) arrival time distribution. The second term captures the physics of multi-phase flow in porous media and is expressed in a relative permeability formulation. Here the overall convective two-phase model is approximated by convolution of the two terms. In order to quantify the effect of the fractional flow curves on production we defined a measure, based on the differences in period of production, that is relatively insensitive to the arrival time. For a number of two dimensional heterogeneous quarter five spot production scenarios and for four fractional flow curves, the two phase behaviour obtained by convolution is compared with the ârealâ two phase flow behaviour (fully 2D numerical solution). Most convolution results are at best in fair agreement with the full physics simulations. This can be attributed to the preferences in the flow direction resulting from total mobility near the fluid front. However, in our results the averaged error in the arrival time prediction of the watercut in the range 5-95% remains smaller than 20% and the error in the oil recovery remains within the 10% range (at 2.7 Pore Volumes injected). In summary, in an extremely fast way, convolution provides us a rough estimate of the full physics displacement.