One outstanding challenge in the physical chemistry of electrolyte solutions is to quantitatively describe the structure and properties of electrical double layers in systems with high surface charges, high ion valencies, high salt concentrations, and spatially varying dielectric permittivity. The classical mean-field Poisson-Boltzmann (PB) theory is physically intuitive and numerically soluble, however, its usage is limited to dilute monovalent salts at low surface charges as it does not account for three essential factors: ion-ion correlations, dielectric variation, and excluded volumes of ions and solvent molecules.A theory addressing these limitations is crucial for understanding many fundamental electrolyte solution phenomena, such as the vapor-liquid interface in ionic fluids, overcharging and charge inversion, repulsion between oppositely charged surfaces, and attraction between like-charged surfaces, among many others. In this thesis, a new electrolyte solution theory is developed that addresses the limitations of the mean-field Poisson-Boltzmann (PB) theory. The validity of this theory is demonstrated by explaining the aforementioned phenomena in a self-consistent and rigorous manner. Going beyond mean-field PB to accurately quantify spatially varying ion-ion correlations, dielectric permittivity, and excluded volume effect is a numerically implausible task. The reason is the need to resolve the electrostatic correlation function at two very different length scales, one associated with ion size (short-range) and the other associated with interface thickness (long-range). Contemporary ways to solve this dual-length scale problem are using a phenomenological approach, a non-local density functional-based approach, or Integral equation-based theories. While phenomenological models have often failed to satisfy well-established Debye-H\"{u}ckel theory in the bulk, the integral equation-based approach, and non-local density functional-based approach use unphysical approximations and non-generalizable weighting functions to reduce the computational cost of this dual-length scale problem. Here, we present a self-consistent field theory entitled, ``Modified Gaussian Renormalized Fluctuation Theory" to overcome the limitations of existing beyond mean-field PB approaches. The main contribution of this work is the introduction of a self-consistent scheme to decompose the correlation function into a short-range contribution associated with the local electrostatic environment and a long-range contribution accounting for the spatially varying ionic strength and dielectric permittivity. This decomposition step makes the dual-length scale problem numerically tractable in a thermodynamically rigorous way. We also account for the excluded volume effect of ions and solvent molecules by including the incompressibility constraint in the partition function. Additionally, we demonstrate the complete numerical method for solving the resultant non-linear equations. We introduce a novel Sturm-Liouville theory-inspired approach that analytically handles the Dirac delta function in the differential equation of electrostatic correlation function, allowing us to employ highly efficient spectral methods. For the problem of vapor-liquid interface in ionic fluids, in the case of symmetric salts, both the coexistence curve and the interfacial tension predicted by our theory are in quantitative agreement with simulation data reported in the literature. We also provide the first theoretical prediction of interfacial structure for asymmetric salt, highlighting the importance of capturing local charge separation. Next, we elucidate the underlying dependence of overcharging and charge inversion on the electrostatic coupling by varying surface charge, counterion valency, salt concentration and dielectric contrast. Consistent with simulations, three characteristic regimes corresponding to weak, moderate, and strong coupling are identified. Important features like the inversion of zeta potential, crowding, and ionic layering at the surface are successfully captured. For weak coupling, there is no overcharging. In the moderate coupling regime, overcharging increases with surface charge. Finally, in the strong coupling regime, ionic crowding and saturation in overcharging are observed. Our theory predicts non-monotonic dependence of charge inversion on multivalent salt concentration as well as the addition of monovalent salt, in quantitative agreement with experiments. We also capture the counter-intuitive phenomena of like-charge attraction and opposite-charge repulsion in multivalent salt solutions and explain their relationship with overcharging. Our theory predicts that the strength of opposite-charge repulsion monotonically increases with salt concentration whereas the strength of like-charge attraction behaves non-monotonically. The addition of monovalent salt to a multivalent salt solution is found to decrease the strength of both opposite-charge repulsion and like-charge attraction. Akin to overcharging, opposite-charge repulsion and like-charge attractions are also outcomes of the heightened ion-ion correlation effect in multivalent ions and there is no inherent causal relationship between overcharging and these two phenomena. Our theoretical predictions for the double layer forces are consistent with the observations reported in experiments and simulations. The thesis concludes by discussing the limitations of our theory and future avenues of research.