Total-internal-reflection deflectometry for measuring small deflections of a fluid surface

We describe a method that uses total internal reflection at the water-air interface inside a large, transparent tank, to measure the interface's deflections. Using this configuration, we obtain an optical set-up where the liquid surface acts as a deformable mirror. The set-up is shown to be extremely sensitive to very small disturbances of the reflecting water surface, which are detected by means of visualising the reflections of a reference pattern. When the water surface is deformed, it reflects a distorted image of the reference pattern, similar to a synthetic Schlieren set-up. The distortions of the pattern are analysed using a suitable image correlation method. The displacement fields thus obtained correlate to the local spatial gradients of the water surface. The gradient fields are integrated in a least-squares sense to obtain a full instantaneous reconstruction of the water surface. This method is particularly useful when a solid object is placed just above water surface, whose presence makes the liquid surface otherwise optically inaccessible.


Introduction
Measuring instantaneous free surface deformations of liquids is of general interest in several practical applications such as in coating and food industries, in large applications such as to study ship wakes, and in off-shore engineering [1,2]. The interest also naturally extends to more fundamental fluid dynamics and physics problems such as studying interfacial fluid instabilities [3,4], droplet dynamics [5,6], wave formation and propagation on the surface of a fluid [7], and in oceanography [8,9].
The methods to quantitatively measure liquid surface behaviour may be broadly divided into two categories based on whether they are intrusive or not. Intrusive methods can be used when the extent of intrusion is small, and the average flow is not significantly disturbed. Traditionally, arrays of resistive (or capacitive) wave probes have been used to study the variation of water level in large setups studying waves [9,10], but can only be installed in sparse distributions separated by gaps of (at least) several centimetres. Less intrusive methods that rely on flow velocities collected using a stereo particle-image-velocimetry setup have also been shown to work for large scale systems [11,12]. Some non-intrusive methods for such measurements, that only use reflections from the water surface and a set of multiple cameras for reconstruction have also been developed [9,13].
A non-intrusive method compatible with smaller, lab scale setups, to resolve deflections of the micrometer to millimeter scale of the free surface, is to use the liquid surface as a refracting or reflecting interface. Usually refraction is used, where the water surface acts as the surface of a lens. A reference pattern is placed underneath the water bath that is contained in a transparent tank. When the light rays from the pattern emerge through the liquid surface, they are refracted due to the jump in refractive index. The variation in heights of the free surface causes further movements of the refracted image of the reference pattern. These movements can be recorded using a camera and analysed to reconstruct the liquid profile. This method is a spin on the well-known Schlieren method, and is known as the free-surface synthetic Schlieren method. It was first proposed by Kurata et al. [14], and since has been matured by the works of Moisy et al. [1] and Wildeman [15] to result in a packaged method that is quick and inexpensive to arrange. The optics of the problem are used to compute the spatial gradients of the liquid surface. The gradient fields are then integrated using a suitable algorithm to obtain a full reconstruction of the imaged area. Even when a fully quantitative reconstruction cannot be obtained, a great deal of qualitative information can be learnt from mere visualisations of the free surface, as used by Fermigier et al. [3], Jain et al. [16] and Chang et al. [5,6]. Non-intrusive acoustic techniques to measure the interface dynamics have also been developed [17].
A few other methods use the reflections from the liquid surface acting this time as a mirror to compute its spatial profile. Cox & Munk [18] were the first to use the specular reflections of the Sun from the sea surface to obtain information about spatial gradients of the water surface. Direct specular reflections can also be obtained from suitably placed lamps, a method used by Rupnik et al. [19] to reconstruct the liquid profile. Another category of such methods uses structured light (such as spatially periodic bright bands of light) that are projected on the free surface. When the surface deforms, the projections also appear distorted. A camera is used to record the movements of the projected fringes, whose phase changes are interpreted to reconstruct the height profile of the liquid surface [20,21]. Such methods have long been used in solid mechanics where extremely small displacements (of the order 10 nm) need to be resolved [22][23][24][25][26]. They have come to be known as 'deflectometry'.
Here we visualise the movements of the water surface by using it as a specularly reflecting surface in a total-internal-reflection (TIR) configuration. Taking inspiration from Moisy et al. [1] and Wildeman [15], we use a fixed pattern, whose distortions by the moving free surface are interpreted in a synthetic-Schlieren sense to obtain displacement fields. Note that contrary to Moisy et al. [1] and Wildeman [15], we use the water surface as a mirror rather than as a lens. From the point of view of a ray-optics problem, the presence of a mirror results in an additional complication as it is the reflecting 'mirror' that undergoes deformation, and not the apparent object that is behind the mirror. We exploit the ray optics in the setup to derive relations between the measured displacement fields and the local spatial gradients of the free surface. Finally we discuss how this gradient information is integrated in a least-squares sense to obtain a fully reconstructed liquid surface profile from the imaged snapshot at a given instant.
The main offering of this particular method is that the liquid surface can be visualised when it is not optically accessible, due to, for instance, the presence of an opaque object above the free surface. An example of such a situation is when a solid projectile is close to slamming onto the liquid surface, and obstructs direct imaging needed for synthetic Schlieren.
As imperfections on a mirror are much easier detected than on a lens, our method is expected and shown to be inherently more sensitive than classical synthetic Schlieren.
The paper is organised as follows: in section 2, we introduce the optics which allow the technique to work, and details of the setup in which we implemented the method. The first stage of the technique involves measuring the displacements of the reference pattern in the mirror plane. The methods to quantify these displacements are discussed in section 3. Next, in section 4, we discuss the relation between these displacements and the deformation of the water surface from which they originate. In section 5, we discuss some subtleties involved in performing the inverse gradient operation in order to finally obtain the final height field, along with an example of the reconstructed surface. In section 6.1 we cover sensitivity, optical limitations and uncertainty estimation. An example of this technique is discussed in section 7, wherein we show a comparison between the measurements and simulations, thereby validating the technique. We end in section 8 with conclusions, the advantages of this technique, and its limitations when compared to other methods which may offer a similar range of accuracy in measurements.

Setup requirements
The setup consists of a water-filled transparent tank with flat walls, a fixed pattern that is allowed to project onto the liquid surface of interest, and a camera to image the reflection from the liquid surface. A light source is used to illuminate the fixed pattern as shown in figure  1.
The light which enters the water tank is refracted towards the interface's normal, as it enters an optically The image from the printed pattern is reflected at the water-air interface and enters a suitably placed high speed imaging camera. At a large enough angle of incidence, the interface acts as a mirror due to total internal reflection, and the camera only captures the mirror image. The light rays illustrate the general optics of the problem. In keeping with the TIR requirement, the angle θ must satisfy: θ ≥ arcsin n a /n w , where n a and n w are the refractive indices of ambient fluid and water respectively. denser medium. Eventually it reaches the air-water interface, where depending on the magnitude of angle of the incidence (represented by θ in figure 1), the light rays might either pass into the surrounding optically rarer medium (here, air) or get specularly reflected as if by a mirror. The latter case is what we aim to obtain, known as total internal reflection (TIR). It requires the angle of incidence at water surface to be greater than the critical angle θ c = arcsin n a /n w , where n a and n w are the refractive indices of air and water respectively. For TIR to occur at an air-water interface, the angle of incidence needs to be greater than θ c ≈ 48.75 • , which may require the water bath depth to be of the order of the lateral width of the tank. Here we use a tank that is 50 cm in length and width, and is filled with water up to a depth of ∼ 30 cm .

Operating conditions
The method described here can be used to visualise the motion of air-water interface only if the light passing from water to air is fully reflected at the surface, which is easily obtained with large incident angles. However, TIR cannot be achieved if the air were replaced by a medium optically denser than water, such as glass (n ≈ 1.52) or silicone oil (n ≈ 1.40): the image of the original pattern (O in figure 1) would always be refracted and never reflected.
With the above conditions satisfied, the air-water surface will only act as a mirror if it exists. Any small contamination floating at the surface disrupts the free surface, such that the 'mirror' disappears at all such locations. This condition also sets the maximum magnitude of deformations that can be measured. Indeed, local and sharp distortions of the air-water interface produce large curvatures. Thus, with the condition θ > θ c still holding true, the light rays reflected at the interface can be deflected away from the sensor of the camera. Additionally, even at small deformations, some raycrossing may occur, especially where curvature is large, making the imaging and interpretation ambiguous.
Note that due to arrangement of the optical setup, the images recorded by an observer at the camera's location are flattened in the y−direction, i.e., along the direction in which light rays are shown to propagate in figure 1 (to the reader, the direction in the plane of the paper). The result is such that a circular object suspended at the water surface appears elliptical. Thus a conversion factor applies to the aspect ratio. This is found by placing a circular disc at the water surface, and measuring the eccentricity of the ellipse that results from the distortion. There is no such distortion along the x-direction (to the reader, normal to the plane of the paper), and the pattern is reflected as is.
As the camera observes the liquid surface laterally, the field of view does not lie in a plane parallel to the camera. To deal with this, sufficient focal depth must be used, and this may demand stronger illumination of the reference pattern. As an alternative, Scheimpflug optics may be used, wherein an additional lens is placed in between the rotated subject-plane and the camera to obtain a corrected properly re-scaled view in the image plane.
Clearly, also other deformations created by optical imperfections in the setup (e.g., curved container walls) can be dealt with using standard digital image correlation techniques performed on the undisturbed image of the pattern.

Quantifying displacement fields
An example of the image of a stationary water surface, as recorded on camera, is shown in figure 2(a). When a disturbance travels across the water surface, it deforms the interface such that the reflected image is distorted, as seen in figure 2(b). The disturbances of the water surface are recorded with time, and the images are processed using an appropriate method to extract displacement vectors u = (u x , u y ), from the movements of the pattern. Here, both u x and u y are functions of the coordinate vectors r = (x, y) that describe the undisturbed liquid surface and, of time t. Two such methods are discussed.

Using cross-correlation
Cross-correlation methods are usually deployed on two subsequent images from a time series (for instance as they are used in particle image velocimetry, PIV), and divide the region of interest into interrogation windows. In typical PIV measurements, a multi-stage algorithm is used, whereby each image is scanned multiple times, with successively decreasing size of the interrogation windows. Cross-correlation techniques, by their very nature, are best used with images that contain a large number of randomly distributed 'particles' (here, dots or squares) [27]. Note that although such a random pattern may be better suited for use with cross-correlation techniques, we here use a pattern with regularly spaced squares due to demanding illumination requirements. Any freely available or commercial PIV program may be used to obtain two-dimensional displacement fields in the x and y directions.
During the interrogation process, we choose window sizes in keeping with the recommendations made by Raffel et al. [27] and Keane & Adrian [28]. However, it can be seen in figure 2(c) that the displacement field can still contain anomalies in some regions. This is due to how the spatial resolution and displacement resolution are affected by the size of the interrogation window. Most of the noise in the data may be smoothened in later stages when reconstructing the water surface (see section 5.2).

Using Fourier Demodulation
When regularly spaced patterns are used (O in figure 1), the images (shown in figure 2) can be processed using Fourier-demodulation (FD) based methods to extract displacement fields. In this case, images from a time series are usually compared to a reference image with the undisturbed pattern. These methods have been commonly used in solid mechanics [23,26] as they can resolve extremely small disturbances which are of use in measuring 2D strain fields. Recently these techniques have been introduced in fluid mechanics [15]. The principle is the following: given a regularly spaced pattern with a periodicity determined by two orthogonal wave vectors k s for s = 1, 2, the intensity profile of the undisturbed pattern, I 0 ( r) is dominated by the Fourier com-ponents corresponding to k s . Here, r is the position vector. A disturbed free surface reflects a distorted pattern, such that the reference intensity profile is slightly deformed, and changes to where u( r) denotes the displacement u of the pattern at position r. By filtering out only the dominant Fourier modes, I 0 (r) transforms into with a s constant. Consequently, the deformed pattern I( r) transforms into for s = 1, 2 , i.e., it is phase-modulated by the disturbances u( r) of the pattern. The latter can be extracted by multiplying g( r) with the complex conjugate of the filtered reference pattern g * 0 ( r) and determining the phase shift For each position r this constitutes a pair of linear equations, which can be readily solved for u( r). An example resulting from this procedure is shown in Figure 2(d). Naturally, some restrictions apply. For example, the components in the signal whose wavelengths are significantly shorter than the pattern wavelength are simply filtered out. The reader can refer to Wildeman [15] for a more detailed discussion on how to select the wave vectors k s of the pattern appropriately.

Comparisons between the two methods
The main difference between using FD and PIV is that while the former compares each image on a stack to a fixed reference image (typically the first in the stack) to calculate the displacement, while the latter involves comparing each image to the preceding one in the series (such that the reference image for a stack is not fixed, and moves along the image stack). Thus when a pattern deforms beyond a certain extent such that no amount of (even distorted) periodicity of the pattern can be detected, the FD method will fail to detect a displacement. In such instances auto-correlation based PIV will still yield a displacement field, which however, will likely contain some inaccuracies.
Since PIV divides the total image into multiple windows, the displacements that occur within the outer margins of the image that are half the width of the interrogation windows, are not resolved. Additionally, the resolution of the displacement field depends on the overlap between adjacent interrogation windows. Obtaining a full-pixel resolution between the image and the displacement field are often computationally very expensive. In contrast, FD yields displacement fields at full-pixel resolution as that of the images being processed, and no information at the margins of the image is lost.
In both methods, displacements may be measured with sub-pixel resolution, but spatial structures smaller than the interrogation window (in the case of PIV), or the wavelength of the pattern (for FD) cannot be easily resolved.

Obtaining surface deformation from projected image distortions
The last task is to relate the displacement vector u( r) to the actual vertical deformation h( r, t) of the liquid sur-face. To do so, we need to consider the ray optics of the setup in some detail. As illustrated in figure 3, an object (source) is placed at position P , from which a light ray travels towards the 'mirror' (here, the air-water interface) from which it is reflected into the camera. Although we measure the displacement fields by tracking the deformation of a fixed pattern, the deformations actually take place at the air-water interface. In other words, it is the mirror that deforms, and makes the image of the object behind it look distorted. The reader is asked to refer to figure 3 as a guide. Since the water surface can either (and often simultaneously) shift over a vertical distance h or tilt by some angle α in either the (y, z)-plane (i.e., parallel to the plane in which camera and pattern lie) or the (x, z)-plane (that is, perpendicular to it), we have here a set of three, generally coupled problems, which we may treat as uncoupled by virtue of the smallness of the free surface deformations that we aim to measure: the 'mirror' may undergo angular Fig. 3 Image reconstructions representing the decoupled 'mirror' deformation problem, which provides the relation between displacements recorded by the camera and the surface deformation h. (a) An angular tilt of the liquid free surface over an angle α in the (y, z)-plane shifts the image P of an object P for the undisturbed free surface to P . A camera looking at the image from the direction in which the lightray (blue) is traveling observes a shift of P indicated by the green arrow. (b) An angular tilt over an angle α in the (x, z)-plane, i.e., perpendicular to the situation depicted in (a), results in a predominantly horizontal shift of the image (green arrow), which again moves from P to P . In this projected view, the camera now is located at the same point as the object P . (c) A vertical shift of the free surface over a distance h in the z-direction results in a vertical shift of the image P over a vertical distance 2h. In the (y, z)-plane of the drawing this leads to a shift observed by the camera as indicated by the green arrow. Clearly, there is no observable shift in the (x, z)-plane. Note that H corresponds to the vertical distance between the reference pattern (the object P ) and the free surface (and thus also to the vertical distance of the image P in the undisturbed mirroring water surface), whereas L denotes the horizontal distance of P to the point where the light ray reflects, i.e., it is related to the angle of reflection θ as tan θ = L/H. The primed quantity H is the distance of P to the deformed free surface (and thus also to the corresponding image P ). The quantities ∆x, ∆y, and ∆z denote the displacements of the image P in space. The vertical displacement h is the quantity that we ultimately want to measure. It is good to realize that all deformations of the free surface are greatly exaggerated for the purpose of this illustration. deflection in two directions (figures 3(a) and 3(b)), or it may simply shift in the vertical direction ( figure 3(c)).
The first step is to relate these three elementary deformations to the shifts that they cause in the camera images, represented by the green arrows in figure 3. The first case, where a tilt over an angle α in the (y, z)plane occurs in isolation, is shown in figure 3(a). A light ray (blue) emerging from P travels towards the mirroring undisturbed free surface (horizontal black line) and reflects towards the camera, the observer. To the observer this light ray appears to travel from a point P , the mirror image of P . With observer fixed, let the mirror tilt by a small angle α. As a result, the image point P now translates to P , which can be found by mirroring P in the tilted free surface (red tilted line). To compute the displacements ∆y and ∆z in the (y, z)plane, we first concentrate on the distance d of the undisturbed and tilted free surface, measured along the line connecting P and P , as indicated in the figure.
On the one hand, d = L tan α, as is obtained from the triangle formed by the endpoints of d and the point where the blue light ray reflects. On the other hand, we can relate the distance H of P to the undisturbed free surface and H of P to the tilted free surface as H = (H − d) cos α = (H − L tan α) cos α, from the triangle formed by H, H and the tilted free surface. Using this relation between H and H, we have where we used that tan θ = L/H and in the approximate equality made use of the fact that α is small such that the trigonometric functions of α can be approximated by their linear Taylor expansions around zero.
Similarly, we have To determine the shift ∆S y,tilt that is observed by the camera (shown by green arrows in figure 3)(a), we need to project ∆y and ∆z onto the plane perpendicular to the viewing direction (i.e., perpendicular to the imaginary light ray represented by the blue dashed line), such that where in the last line we have approximated α by the local slope of the free surface at the point where the light ray touches the interface: ∂h/∂y = − tan α ≈ −α + O(α 2 ).
The second case, depicted in figure 3(b), corresponds to a tilt over an angle α in the (x, z)-plane, where for clarity the light ray has not been drawn since in this projection it would coincide with the line connecting the object P with the free surface. This case may be analysed in a very similar manner as the first. The displacements ∆x and ∆z in the (x, z)-plane can now be directly deduced from the orthogonal triangle formed by P , P and the intersection of the vertical through P and the horizontal through P , together with the relation H = H cos α obtained from the lower triangle in the figure ∆x = 2H sin α = 2H cos α sin α ≈ 2Hα + O(α 3 ) , (8) and where it is good to note that the latter is of order α 2 and also does not lead to a shift in the image plane of the camera. Therefore, the shift ∆S x,tilt that is observed by the camera equals again approximating α by the local slope of the free surface: ∂h/∂x = − tan α ≈ −α + O(α 2 ). The third case corresponds to a vertical shift h of the free surface in the positive z-direction as depicted in figure 3(c). Clearly such a shift only leads to a corresponding shift of the camera image in the (y, z)-plane, where the displacement of the image P to P is also a simple vertical shift over a distance 2h, that is ∆z = 2h. As a consequence the shift ∆S y,shift that is observed by the camera, i.e., the projection of ∆z onto the plane perpendicular to the viewing direction is equal to where we took into account the opposite direction of the shift as compared to that of the first case by the minus sign. Finally we note that by using the same symbol (h) for the vertical shift and the local vertical deformation of the free surface, the result is already stated in terms of h(x, y, t). In summary, we find that our elementary deformations of the free surface result in a displacement field (∆S x , ∆S y ) in the camera image that is given by: where h(x, y, t) is the deformation of the free surface that we are after.
The second step is to relate the displacement field ∆ S = (∆S x , ∆S y ) in the camera image to the field u discussed in the previous section, which is a rather subtle one. One may be tempted to just equate the two, but then one overlooks that structures on the free surface appear deformed in the camera image since the latter is observing the free surface under an angle θ. E.g., circles on the free surface appear like ellipses with their short axis in the y-direction in the camera image. Naturally, one will correct the camera images for these kind of deformations, but now one has two options, namely to perform the FD or PIV analysis either before or after this correction, i.e., one may determine the displacement field u either before or after correcting the camera images. In general we have found it advantageous to first correct the grid in the camera images, such that the coordinates r = (x, y) in the corrected camera images correspond to the coordinate system to the free surface (conveniently denoted by the same symbols and corresponding to the notation that has been used throughout the article).
For the case depicted in figure 3, transforming back to the coordinate system attached to the free surface amounts to dividing the y-coordinate of the camera image by cos θ. Since the displacement field (∆S x , ∆S y ) has been related to the vertical deformation field h(x, y, t), and u(x, y, t) has been obtained in the coordinate system attached to the free surface, we need to make the same transformation for the y-component of ∆ S, i.e.
In the case that the field of interest of the free surface is small compared to the distance H (and subsequently also to L), it may well be sufficient to assume that θ is constant and that no correction in the x-direction is necessary, as in equation 13. In general however, θ is not constant, but a function of x and y, and a similar correcting factor may also be necessary in the x-direction for those points that are far away from the center.
In any case, one may rewrite equation 13 as whereĵ denotes the unit vector in the y-direction and ũ = (ũ x ,ũ y ) and T are defined as (ũ x ,ũ y ) = u x , u y cos 2 θ and T = cos 2 θ tanθ = 1 2 sin(2θ) .
This equation constitutes a partial differential equation for the vertical deformation field h( r, t) in terms of the experimentally known (ũ x ,ũ y ) and will be the basis of our analysis in the coming sections. It is convenient to split equation 14 into components: On careful observation of figure 3 one may notice that, drawing a light ray from the displaced image point P towards the camera also causes a shift of the point where it reflects from the free surface. That is, one is not exactly measuring the vertical shift and angular tilt of the free surface in the point (x, y) but in a slightly shifted point. In principle one could correct for such an image shift, but if h and especially ∇h are not varying too quickly on the free surface one can neglect this effect. Since ray crossing limits the second derivatives of h with respect to x and y, as will be discussed in detail in section 6.1, this condition is generally fulfilled.

Recasting the integrand using an integrating factor
Note that equation (14) cannot be directly integrated due to the additional dependence on h. Thus we recast the expression using an integrating factor, under the assumption that T is constant over the region of interest, i.e., independent of x and y. Equation (17) can be re-written as Similarly, equation (16) can be re-written using the same integrating factor Equations (18) and (19) can be combined using vector notation as or, The gradient fields in x and y directions, that are to be integrated over, are expressed in the form shown on the right hand side of equation (21). The result obtained from surface integration is divided by the factor exp( yT H ) to obtain the final height field h(x, y). With equation 21, we have now recast our original problem in a conservative form where ξ is the known vector field, and f is to be determined. Mathematically such an expression can be directly integrated since ∇ × ξ = ∇ × ∇f ≡ 0. However, since ξ is only approximately known due to unavoidable noise in the experiments, some additional care is needed during the integration.

Inverse gradient operation
The inverse gradient operation is performed on equation (21) to obtain the final result where f 0 is an integration constant, connected to the absolute height of the free surface. In the following discussion, f 0 is set to zero for convenience. One way to integrate over the gradient information ξ is to start at a reference point (x r , y r ), and integrate along a path such that However, using this approach, any noise in the local gradient information may get added over the path of integration [1]. Moreover, in a discretised implementation of this method, it is not clear how the final result would be modified if the order of integration along the paths in x and y direction were switched. Both drawbacks can be avoided by using a 'global' approach. This is done by building a linear system of equations, replacing the gradient by a 2nd-order centred finite difference operator. If the (x, y) space is discretised by M ×N elements, there are M N variables to be determined (corresponding to the discretised height field f (x, y, t)), while there are 2M N knowns (the gradient information stored in u) in the system, leading to as many equations. Thus, we are dealing with an over-determined matrix system, which cannot be simply inverted. The inversion is therefore performed while minimising a residual cost function. More details can be found in [1,29].
The least-squares solution thus found has the effect of smoothening out local outliers present in the gradient fields. An efficient MATLAB implementation was written and made public by D'Errico [30]. More details on global least squares reconstruction, and further advanced methods can be found in the works by Harker & O'Leary [29,31,32]. We use the implementation by D'Errico which is now commonly used in reconstruction problems that involve an inverse gradient operation to be performed on a mesh of spatial gradients [1,[33][34][35]. An example of the reconstructed surface profile, based on the typical displacement field shown in figure 2(d), is shown in figure 4. A more systematic experiment, along with comparisons with simulations is discussed in section 7.
Now, finally, one may ask what can be done when T is not constant, that is, when the region of interest at the free surface is not small compared to H. In principle one may use an integration factor where the exponent is an integral over T /H, which would however add additional complexity to the analysis. One may however also ask how large of an error one makes by approximating T /H by a constant. It turns out that this error is relatively small, since the relative sensitivity of the method to a vertical shift is small, as will be discussed in detail in the next section. Fig. 4 Reconstructed surface profile of water from the displacement field shown in figure 2(d). The arbitrary disturbances on the water surface were recorded and measured over a small section of the total water surface in the bath, that is shown above.
6 Sensitivity, limitations and error estimation

Sensitivity
Starting from equation (14) which relates the surface profile height h(x, y, t) to the measured displacement field ũ(x, y, t), one immediately realizes that there are two manners in which the surface profile height may influence the displacement field, namely by a tilting of the interface (corresponding to the first term on the right hand side, ∼ ∇h) or by a vertical shift (the second term which is proportional to h). We will now address the sensitivity of the setup, where we will start with assessing the relative sensitivity of a tilt versus a vertical shift.
Since tilt and shift are usually correlated, we start by performing a modal decomposition of the surface height profile, where it suffices for our purposes to concentrate on the y-direction only Rewriting the y-component (17) of equation (14) as and inserting equation (25) yields, for each of the modes separately −2HkA k cos(ky)−2T A k sin(ky) =ũ y,tilt +ũ y,shift . (27) Clearly, the two terms of these equation do not attain their maxima in the same points, as a result of the fact that sin(ky) is zero where its derivative is maximal and vice versa, but one may easily compute the respective maxima and determine the relative sensitivity as the ratio of thesẽ Note that this ratio is independent of the amplitude A k . Since the wavelength of even the largest structures that are to be observed is usually much smaller as the distance of the liquid surface and the pattern, i.e., λ H, the above ratio is typically much smaller than one, which implies that the setup is much more sensitive for a tilting of the surface than for a vertical shift u y,tilt u y,shift . To put this difference in absolute terms, we note that the detection of the displacement field u is bounded by the sensitivity of the method used to obtain it which provides a minimum detectable displacement δv y,min which is some fraction of the pixel size of the measured image. Usingũ y,tilt > δũ y,min , we find that 2HkA k δũ y,min , yielding A k δũ y,min λ 4πH .
From the above we can immediately conclude that the deformations that are visible with our method are much smaller than the spatial resolution of the displacement pattern. For the example of Fig. 4, where H = 30 cm, and the typical wavelength of the structure is λ ≈ 2 cm, we find that λ/(4πH) ≈ 0.005. Using a spatial resolution δu y,min = 100 µm, we obtain that the minimal displacement δh min,tilt that is discernible through the detection of the tilted interface equals δh min,tilt = 0.5 µm, and that this sensitivity may (at least theoretically) be increased by increasing the distance H between camera/pattern to the liquid surface. Similarly, we obtain for the sensitivity for a vertical shift that 2 tan θA k δũ y,min , or, As expected, the result is independent of the wavelength and much larger than it is in the case of a tilted interface. In fact, using the same spatial resolution in the case of the example of Fig. 4 (θ ≈ 55 • ) we have δh min,shift = 50 µm, i.e., the setup is two orders of magnitude less sensitive for a vertical shift than for a tilt.
Conversely, this means that if two patterns differ by a vertical shift, i.e., h 1 (x, y, t) = h 2 (x, y, t) + ∆h(t), the difference between h 1 (x, y, t) and h 2 (x, y, t) would be very difficult to detect, especially if ∆h is of the same order as h 1,2 , which would usually be the case in experiment. Here, the contribution of ∆h to the signal would be typically two orders of magnitude smaller than that of the surface deformation features. This implies that, even in a time series, there may be a shift between the profiles determined at different moments in time that is extremely hard to detect, if at all. This makes the method most suitable in the case that there exists a reference point on the interface where no deformation is expected.

Limitations
The setup has several limitations originating from the fact that it makes use of the liquid surface as a deformed mirror, which we will discuss in sequence in this subsection.

Mirroring condition
Total internal reflection will only happen if the angle of incidence φ i on the deformed liquid surface is larger than the minimal angle φ i,min for which total internal reflection will take place, i.e., φ i > φ i,min = arcsin n a n l .
Now the angle of incidence is determined by the angle θ at which we look at the pattern and the slope of the liquid surface in the y-direction ∂h/∂y, namely φ i = θ − arctan(∂h/∂y) which limits the slope to ∂h ∂y tan θ − arcsin n a n l , or, A k λ 2π tan θ − arcsin n a n l .
As long as the typical length scale on which the pattern changes (λ) is sufficiently larger than the amplitude (A k ) we seek to measure, satisfying the above condition will not be a serious problem, provided θ is not chosen too close to arcsin(n a /n l ).

Ray crossing
Two incident, parallel rays will cross before reaching the camera if the local radius of curvature R of the liquid surface is smaller than the distance of the camera H/ cos θ. Since for small deformations the radius of curvature can be approximated as 1/R ≈ ∂ 2 h/∂y 2 , we obtain, using modal decomposition (25) For the example of Fig. 4 (H = 30 cm, λ ≈ 2 cm, θ ≈ 55 • ), this will lead to A k 19.5 µm. This is quite a stringent requirement, which can be improved by moving the camera closer to the liquid surface, or decreasing H. As discussed above, doing so will however lead to a loss of sensitivity. From another perspective, the necessary condition to prevent ray crossing, 1/R cos θ/H, sets an upper limit to the second order spatial derivatives of h, which implies that h should vary little on the length scale set by h itself. This implies that a shift of the result in the order of the measured amplitude, as discussed at the end of section 4, will negligibly impact the reconstructed free surface deflection h.

Error estimation
The method is prone to some systematic and random errors that in the end will propagate into the measurement result, the deformation of the interface h(x, y, t). Some of those are quite generic for systems making use of high-speed optical image acquisition, and find their origin in the specifications of the camera (spatial and temporal resolution, motion blur, pixel sensitivity) and have to be addressed by using a camera that is suitable for the particular problem at hand [36]. Others are related to the quite elaborate image data processing to first detect the displacement field u(x, y, t) in the image plane (using PIV or FD) and to secondly compute h(x, y, t) with the spatial integration method, and are difficult to assess or control. Here it is crucial to employ a scheme that integrate the displacement field in a global least square sense (as discussed in Subsection 5.2), as otherwise especially systematic errors may be cumulatively integrated and lead to substantial errors in h.
Relevant from the perspective of the current setup is how errors in the main parameters H and θ propagate in the final interface profile h. Based upon the sensitivity results of Subsection 6.1 one may expect that the influence of errors in H are more significant than those in θ. More quantitatively, we may use the modal decomposition (25) in equation (13) to determine how a variation ∆H in H propagates into a variation ∆A k of the amplitude A k of mode k, leading to where we have used that the wavelength of the observable structures are usually much smaller than H (i.e., λ/(2πH) 1, such that the second term in the denominator is small everywhere except close to where the slope of the interface is zero. Similarly, we can write for the propagation of a variation ∆θ in θ that where the dominant term (for cot(ky) not too small) has been kept in the second approximation. The first term is much smaller than one whereas the second is typically of order unity, such that the relative error in θ is multiplied by a small number. This is good to realize when setting up the experiment: it is more crucial to assure that the pattern is positioned such that H can be considered constant over the region of interest, and some compromise in the constancy of the value of θ can be made in order to reach that goal.
7 Example and validation: Water surface deflection due to air cushioning under an approaching plate Validation of the experimental method is difficult due to the sensitivity of the method. When one tries to use known or macroscopically observable menisci around immersed objects, the problem is that the interface disturbance close to the object is not observable due to the large local deflection and curvature. This implies that one may only observe the far-field exponential decay which is hard to relate to a physical length scale. This leaves the observation of water waves (as has been done qualitatively in earlier Sections) or the deformation of the interface due to the impact of an object. We will now turn to the latter and, for the purpose of validation reproduce some results from [38] in figure 5. The experimental setup is described in figure 6(a): a flat disc is slammed onto a stationary water bath with a controlled velocity. The approaching disc pushes out the ambient air from the gap in between itself and the water surface. The stagnation pressure set up under the disc centre deflects the water surface away. The (azimuthally averaged) measured profiles are shown in panel (b) at various times before impact (τ ). The measurements at r = 0 are compared with two-fluid boundary integral simulations described in [37,[39][40][41][42]. The favourable comparison indicates that the measurement technique is successful at resolving deflections of the order of micrometres up to several tenths of millimetres. For additional information, we refer to [38]. , non-dimensionalised using inertial scales D and V , and compared with the boundary integral simulation from Peters et al. [37] and [38].

Conclusions
We described a TIR-based method to measure smallscale deformations of a water surface, consisting of two steps: First, the movement of the water surface is measured by recording the deformation of a reference pattern that is reflected in the water surface. The displacement of the reference pattern is then quantified using an image correlation method such as PIV or FD.
Secondly, these displacements are interpreted as projections in the two-dimensional image plane, and related to the instantaneously deforming water surface and its spatial gradients. By decoupling the light paths when the reflecting surface either undergoes an angular deflection, or a vertical translation, we build a system of equations that relate the pattern deformation to the local surface deflection. This second step thus involves recasting the measured displacement fields to a suitable integrable form, and calculating the final height field.
Since the image manipulation and subsequent solution may become quite complex, it is wise to test the setup using an axisymmetric deformation of the free surface, before turning to the measurement and analysis of less symmetric situations. Some images, the displacement fields obtained from them, and an example code for reconstruction are provided in supplementary material.
A relative drawback of TIR-deflectometry arises from the high sensitivity it offers: it requires the water surface to be very well isolated from external sources of noise. This high degree of isolation from mechanical disturbances limits the method's application to wellcontrolled environments. Another consequence of the sensitivity is that using menisci of a stationary object for calibration purposes is difficult, since deflections easily become too large to be measurable.
An application of this method was discussed by measuring the water surface deflections due to air-cushioning under a plate that is about to slam on it. Good comparison of the measurements with boundary integral simulations validate the technique for measurements up to tens of micrometres. Some more examples of the use of this method are described in ref. [43, chapter 6] by measuring micron-scale waves on a water surface, and showing successful comparisons with a theoretical model, thus showing its effectiveness in resolving precise micron scale deformations.
The method's greatest merit lies in it using total internal reflection at the water surface. This implies that whatever moves above the water surface remains invisible to the camera. Additionally, sub-micron resolution of the interface deflections is readily achieved.