Solvation is the interaction of solute and solvent. Every biological interaction hap- pens in a solvent. Most technical procedures occur within solvents and geological processes are too mediated by their solvents. Understanding these effects of solva- tion is therefore critical for the understanding of biology, technology, and geology. While static properties of solvation shells, like coordination numbers and radial dis- tribution functions, are well understood, the dynamic properties of these open sys- tems are rarely studied. Furthermore, an interesting solvation-based phenomenon is the separation of the geochemical twins Zirconium and Hafnium in fluoride-bearing media. In this work, I present a method for evaluating Markov models of solvation shells and investigate ways of combining the solvent models with solute models. Moreover, I attempt to find a suitable two-site and three-site model of hydrogen flu- oride, specialized for its interaction with metal ions. By simulating aqueous and pure HF for several combinations of q, σF , and ϵF and evaluating density, peaks of radial distribution functions as well the solvation free energy of NaF in HF, I hoped to find a suitable model. A three-site model for HF is parameterized by recreating the electrostatic potential of HF with a classical force field, focusing on the location of maximum potential which takes a conical shape around the tip of the ellipsoid and is not located at the poles. A new method is presented which allows the automatic detection of coordination polyhedra based on reference structures and Steinhardt-order parameters. The Lennard-Jones parameter space for tetravalent cations is explored and analyzed in terms of static solvation shell quantities. Finally, the thermal contraction of the solvation shells of Zr4+ and Hf4+ in 1 M HF was investigated using classical MD simulations. The Markov models of solvation shells indicated that solvent dynamics couple close to the solute. Additive combined models yielded slightly higher timescales compared to their individual components. The opposite is true for the multiplicative models which performed just as well or even worse than their components. The parameterization of HF, for the two-site model, yielded two parameter combi- nations that could reproduce three of the five target quantities, the relevant peaks of the F-H and H-H radial distribution functions, as well as the solvation free energy of NaF in HF. After choosing a topology for the three-site model, the Lennard-Jones parameter scans were unable to yield stable simulations of aqueous HF. The project was therefore discontinued and I settled for a recently published HF model. The parameterization of metal cations yielded a very robust result. Static solvation shell properties exist on continuous regions of similar value in the parameter space. These regions appear as a diagonal lines in the log(ϵ M ) − σM parameter space. This behavior is also observed in the coordination polyhedra found by the novel method. The thermal contractions of solvation shells of tetravalent cations could be observed for the four ionic ligands. The contractions are a result of water molecules increas- ing their distance to the central cation. Their missing repulsive Coulomb interaction allows the ionic ligands to move in closer to the central cation, thus causing the thermal contraction. Furthermore, we observed a two-state system for the solvation shells at high temperatures which consists of octahedral and tetrahedral solvation shells interchanging each other. The herein presented results offer new methods for analyzing solvation shells. Firstly by constructing Markov models of these open systems to study their dynamics. Secondly, by automatically determining the coor- dination polyhedron, which is essentially an analysis of the angular distribution of solvent molecules in the solvation shell. The parameterization attempts of HF depict the difficulty of finding parameter com- binations that match all fitting targets, albeit the searched parameter space was rather small. The parameter space for tetravalent cations shows an extremely robust result which can yield the basis for future parameterization attempts. Finally, the peculiar thermal contractions could be explained through classical MD simulations. This shows the power of this method for studying hard ionic systems.