Elucidating multiscale, multiphase and multiphysics phenomena of flow and transport processes in porous media is the cornerstone of numerous environmental and engineering applications. Several factors including spatial and temporal heterogeneity on a continuity of scales, the strong coupling of processes at such different scales at least at a localized region within the domain, combined with the nonlinearity of processes calls for a new modeling paradigm called multiscale models, which are able to properly address all such issues while presenting an accurate descriptive model for processes occurring at field scale applications. Furthermore, the typical temporal resolution used in modern simulations significantly exceeds characteristic time scales at which the system is driven and a solution is sought. This is especially so when systems are simulated over time scales that are much longer than the typical temporal scales of forcing factors. In addition to spatial and temporal heterogeneity, mixing and spreading of contaminants in the subsurface is remarkably influenced by oscillatory forcing factors. While the pore-scale models are able to handle the experimentally-observed phenomena, they are not always the best choice due to the high computational burden. Although handling across-scale coupling in environments with several simultaneous physical mechanisms such as advection, diffusion, reaction, and fluctuating boundary forcing factors complicates the theoretical and numerical modeling capabilities at high resolutions, multiscale models come to rescue. To this end, we investigate the impact of space-time upscaling on reactive transport in porous media driven by time-dependent boundary conditions whose characteristic time scale is much smaller than that at which transport is studied or observed at the macroscopic level. We first introduce the concept of spatiotemporal upscaling in the context of homogenization by multiple-scale expansions, and demonstrate the impact of time-dependent forcings and boundary conditions on macroscopic reactive transport. Proposing such a framework, we scrutinize the behavior of porous media for ``quasisteady stage time'' (thousands of years), where there is an interplay between signal frequency and the three physical underlying mechanisms; advection, molecular diffusion and heterogeneous reaction. To this end, we demonstrate that the transient forcing factors augment the solute mixing as they are combined with diffusion at the pore-scale. We then derive the macroscopic equation as well as the corresponding applicability criteria based on the order of magnitude of the dimensionless Peclet and Damkohler numbers. Also, we demonstrate that the dynamics at the continuum scale is strongly influenced by the interplay between signal frequency at the boundary and transport processes at the pore level. We validate such a framework for reactive transport in a planar fracture in which the single-component solute particle is undergoing nonlinear first-order heterogeneous reaction at the solid-liquid interface, while the medium is episodically influenced by time-dependent boundary conditions at the inlet. We also present the alternative effective transport model at a much lower cost, albeit at the regions where the corresponding applicability criteria are satisfied. We perform direct numerical simulations to study several test cases with different controlling parameters i.e. Peclet and Damkohler numbers and the space/time scale separation parameters; the ratio of characteristic transversal and longitudinal lengths $\varepsilon$ and $\omega$; the ratio of period of time-fluctuating boundary conditions to the observation time scale. A rigorous justification of the effective transport model for the given applicability conditions is demonstrated, essentially by comparing the local vertically averaged microscopic simulations with their corresponding macroscopic counterparts. Moreover, as a special case, we employ a singular perturbation technique to look at the effective model for vertical mixing through a narrow and long two-dimensional pore. We obtain explicit expressions for dispersion tensor as well as the other effective coefficients in the coarse-scale homogenized equation. Our analysis manifests robustness of the sufficient and necessary applicability constraints which validate the upscaled model as a solid replacement of the pore-scale one within the accuracy prescribed by homogenization theory. While a deterministic model is sufficiently robust for a plethora of subsurface applications, a more realistic setting is often required when dealing with other scopes of engineering applications, e.g. reservoir engineering and enhanced oil recovery. Rigorous modeling of these systems calls for sophisticated strategies for uncertainty quantification and stochastic treatment of the system under study. Such an uncertainty is inherent to, and critical for any physical modeling, essentially due to the incomplete knowledge of state of the world, noisy observations, and limitations in systematically recasting physical processes in a suitable mathematical framework. To this end, accurate predictions of outputs (e.g. saturation fields) from reservoir simulations guarantee precise oil recovery forecasts. These quantitative predictions rely on the quality of the input measurements/data, such as the reservoir permeability and porosity fields as well as forcings, such as initial and boundary conditions. However, the available information about a particular geologic formation, e.g. from well logs and seismic data of an outcrop, is usually sparse and inaccurate compared to the size of the natural system and the complexity arising from multiscale heterogeneity of the underlying system. Eventually, the uncertainty in the flow prediction can have a huge impact on the oil recovery. Consequently, we also develop a probabilistic approach to map the parametric uncertainty to the output state uncertainty in first-order hyperbolic conservation laws. We analyze this problem for nonlinear immiscible two-phase transport (Buckley-Leverett displacement) in heterogeneous porous media in the presence of a stochastic velocity field, where the uncertainty in the velocity field can arise from the incomplete description of either porosity field, injection flux, or both. Such uncertainty leads to the spatiotemporal uncertainty in the outputs of the problem. Given information about the spatial/temporal statistics of the correlated heterogeneity, we leverage method of distributions (MD) to derive deterministic equations that govern the evolution of single-point CDF of saturation in the form of linear hyperbolic conservation laws. We first derive the semi-analytical solution of the raw CDF of saturation at a given point, for the cases in which two shocks are present due to the gravitational forces. Then, we describe development of the partial differential equation that governs the evolution of the raw CDF of saturation, subject to uniquely specified boundary conditions in the phase space, wherein no closure approximations are required. Hereby, we give routes to circumventing the computational cost of Monte Carlo scheme while obtaining the full statistical description of saturation. This derivation is followed by conducting a set of numerical experiments for horizontal reservoirs and more complex scenarios in which gravity segregation takes place. We then compare the CDFs as well as the first two moments of saturation computed with the method of distributions, against those obtained using the statistical moment equations (SME) approach and kernel density estimation post-processing of exhaustive high-resolution Monte Carlo simulations (MCS). This comparison demonstrates that the CDF equations remain accurate over a wide range of statistical properties, i.e. standard deviation and correlation length of the underlying random fields, while the corresponding low-order statistical moment equations significantly deviate from Monte Carlo results, unless for very small values of standard deviation and correlation length.