To facilitate computational investigation of intermolecular interactions in the solution phase, we report the development of ALMO-EDA(solv), a scheme that allows the application of continuum solvent models within the framework of energy decomposition analysis (EDA) based on absolutely localized molecular orbitals (ALMOs). In this scheme, all the quantum mechanical states involved in the variational EDA procedure are computed with the presence of solvent environment so that solvation effects are incorporated in the evaluation of all its energy components. After validation on several model complexes, we employ ALMO-EDA(solv) to investigate substituent effects on two classes of complexes that are related to molecular CO₂ reduction catalysis. For [FeTPP(CO₂-κC)]^2- (TPP = tetraphenylporphyrin), we reveal that two ortho substituents which yield most favorable CO₂ binding, -N(CH₃)₃ ⁺ (TMA) and -OH, stabilize the complex via through-structure and through-space mechanisms, respectively. The coulombic interaction between the positively charged TMA group and activated CO₂ is found to be largely attenuated by the polar solvent. Furthermore, we also provide computational support for the design strategy of utilizing bulky, flexible ligands to stabilize activated CO₂ via long-range Coulomb interactions, which creates biomimetic solvent-inaccessible "pockets" in that electrostatics is unscreened. For the reactant and product complexes associated with the electron transfer from the p-terphenyl radical anion to CO₂, we demonstrate that the double terminal substitution of p-terphenyl by electron-withdrawing groups considerably strengthens the binding in the product state while moderately weakens that in the reactant state, which are both dominated by the substituent tuning of the electrostatics component. These applications illustrate that this new extension of ALMO-EDA provides a valuable means to unravel the nature of intermolecular interactions and quantify their impacts on chemical reactivity in solution.