Coastal aquifers serve as the major source of freshwater for millions of people living in coastal areas. Seawater Intrusion (SWI), the advancement of saltwater into coastal aquifers due to convective flow caused by density difference between freshwater and saltwater, can lead to salinization of freshwater resources and adversely affecting the lives of the residents of coastal zones. A great number of coastal aquifers worldwide contain significant fracture networks. These fracture networks can act as preferential pathways for fluid flow and transport of dissolved chemicals such as salt leading to the intensification of SWI. Therefore, understanding the effects of fracture characteristics and geometry on the extents of SWI is of high importance for researchers working on different realms of geosciences. Numerical simulation has been proven to be an efficient tool to investigate SWI in fractured coastal aquifers (FCAs). However, due to the complex nonlinear nature of the problem, the inter-relation of associated parameters and the model output is not completely clarified in the literature.

The aim of this study is to use polynomial chaos expansion (PCE) as a tool for performing uncertainty analysis on SWI in FCAs simulated using the coupled discrete fracture network (DFN) and variable-density flow (VDF) models. Applying DFN-VDF requires detailed description of fracture parameters, most importantly the geometric characteristics and hydraulic properties of the major fractures. The uncertainty associated with the identification of such parameters can result in major effects on the model outputs and thus the validity and capability of the model. In order to investigate and quantify the effects of these uncertainties on the model outputs, global sensitivity analysis is performed on different sets of simulations with different fracture characteristics. Two fractured configurations of Henry Problem are introduced and investigated throughout the process: (i) Single horizontal or vertical fracture configuration and (ii) Network of orthogonal fractures configuration. The DFN-VDF models are simulated and solved using finite element based framework of COMSOL Multiphysics®. Prior to the application of GSA on the simulation sets, the validity and robustness of COMSOL on handling coupled DFN-VDF models was checked and compared against an in-house research code where fractures are modeled by considering heterogeneity of material without reduction of the dimensionality. Then, for the mentioned scenarios, GSA is performed in 3 steps: (1) experimental design set is generated using quasi Monte-Carlo sampling technic, (2) then, to increase the computational efficiency of the approach, sparse PCEs are built for each model output, and finally (3) Sobol' indices, which are directly derived from the surrogate PCEs coefficients, are used as sensitivity measures to investigate the primary sources of uncertainties in the model outputs.

The results show that the salinity distribution is sensitive to the fracture hydraulic conductivity. It was observed that for single horizontal fracture configuration, the location of the fracture also plays a significant role in generating uncertainties on all designed metrics. The imperfect knowledge of fracture location and spacing affects mainly the toe position and total flux of saltwater entering the domain of the aquifer. An important aspect of this work is contributed to calculation and demonstration of marginal effects based on the obtained PCEs to understand the effect of fracture characteristics on SWI. The results can be used as a tool for technical and managerial purposes associated with monitoring, control and prevention of SWI in FCAs.