Numerical modeling of groundwater flow and mass transport is increasingly becoming a reference criterion nowadays for water resources assessment and environmental protection. To render the model reliable for future predictions, the model structure and parameters have to be characterized as close to the reality as possible. The process of model identification by integrating measured parameters and observed model states is so called inverse problem. A series of methods has been proposed to solve the inverse problem in the past several decades and its evolution is discussed in the thesis. The main point of this thesis is to propose two stochastic inverse methods to estimate model parameters which cannot be described with a Gaussian distribution, e.g., hydraulic conductivities, integrating nonlinear model observations, e.g., hydraulic head data. The first is the normal-score ensemble Kalman filter (NS-EnKF) constructed on the basis of the standard EnKF. The standard EnKF is widely used as a real time data assimilation technique due to its advantages, e.g., computation efficiency and ability to assess model uncertainty. However, it is known to perform optimally when the model parameters and state variables follow multiGaussian distributions. To extend the application of the EnKF to nonGaussian distributed state vectors, such as those in channelized fluvial-deposit aquifers, we introduce the normal-score transformation into the EnKF, to propose the NS-EnKF. The augmented state vector consisting of model parameters and state variables are normal-score transformed so that they follow marginal Gaussian distribution. Then, the transformed vectors serve as input to the EnKF, which now operates on marginally Gaussian distributed variables. The updated vectors are then back transformed to the original distribution space. The effectiveness of the proposed method is assessed in two synthetic bimodal aquifers, where the NS-EnKF is found to perform better than the standard EnKF in characterizing the bimodal structure of the hydraulic conductivities and in the subsequent flow and transport predictions. The second method is a pattern search-based inverse method (PSINV), a novel inverse method that goes beyond the minimization of an objective function. The unknown model parameters are simulated by searching for pattern matches including parameters and state variables arranged in pre-established spatial templates. The attractive features of the proposed method, PSINV, include that (i) the estimated parameters do not necessarily follow a normal distribution, avoiding the pitfalls associated with the multiGaussian distribution, such as the lack of connectivity on extreme values; (ii) the model structures and the relationship between model parameter and state observations are characterized by multiple-point geostatistics, which renders it possible to apply the method in highly heterogeneous aquifers of complex geometries. The performance of the method is evaluated in two synthetic reservoirs composed of two facies, sand and shale.