One of the easiest ways to interpolate a new value between two known values is linear interpolation. It takes the closest values and creates a linear function between them.
In one dimension, linear interpolation is pretty trivial. For example, let's assume that we have 20 random points $P$ between 0 and 1 (\textcolor{Red}{\textbullet} and \textcolor{Blue}{\textbullet} in Figure \ref{fig:one-dim-interpolation}) and have a new point $I$ (\textcolor{Green}{\textbullet}) at $0.4$ for which we want to interpolate. Finding the two closest points (\textcolor{Red}{\textbullet}) above and below is trivial as there is only one dimension to compare. Now, if we have measured a value $f(P)$ for each of these points, a straight line (\textcolor{LightGreen}{\textbf{|}}) between the two closest values can be drawn and an interpolated value for $f(I)$ can be found.
In two dimensions things get more complicated as we now have a set of points with $X$ and $Y$ coordinates (Figure \ref{fig:3dinterpolate-1}). One fast way to find the closest points to the point that should be interpolated is using Delaunay triangulation. This is a method to separate the space between the points into triangles while trying to maximize their smallest angle and to make sure no other point is inside the circumcircle of the triangles\footcite{Delaunay}. Connecting the centers of the circumcircles results in a Voronoi diagram.
Afterwards, the closest three points can be found very quickly by checking the nodes of the surrounding triangle (Figure \ref{fig:3dinterpolate-2}). If we now again have a function $f(X,Y)$ similar to the one-dimensional example, we can create a unique plain through the three points and get the interpolated value for any pair of $X$ and $Y$ on this layer. (Figure \ref{fig:3dinterpolate-3})
This approach has the advantage that it can be extended in more than two dimensions by replacing the triangle in the Delaunay triangulation with an n-simplex in n dimensions. The \texttt{scipy.spatial.Delaunay} python function allows to quickly calculate it thanks to the \texttt{Qhull} library\footnote{\url{http://www.qhull.org/}}. One noticeable limitation of this method is that data can't be extrapolated. Therefore, the possible output is limited to the convex hull of the input parameter space (as seen in Figure \ref{fig:3dinterpolate-2}).
For doing the actual interpolations, the \texttt{scipy.interpolate.griddata} function is used with the \texttt{method="linear"} argument which itself uses \texttt{scipy.interpolate.LinearNDInterpolator} to do an interpolation similar to the one described above. The function is given a $6\times n$ matrix of the six parameters and an $n$-sized list of the water retention fraction for those $n$ simulations. In addition, \texttt{griddata} supports not only calculating interpolations for one set of parameters, but also for lists of parameters which allows to quickly generate 2d diagrams as seen in Figure \ref{fig:griddataresults}.
Most notable about the results of the griddata interpolation (Figure \ref{fig:griddataresults}) are the many fine details that can be seen. This is mostly caused by the fact that this method only uses the closest values for interpolations and therefore there is no smoothing. These details might just be random derivations of the simulation and not a higher resolution of the data. Another thing that can be seen in the bottom right corner of Figure \ref{fig:griddata1} is that griddata can't extrapolate data.
