Digital Image Correlation for stereoscopic measurement of water surface dynamics in a partially filled pipe

This paper presents a novel method for measuring three-dimensional (3D) water surface dynamics in a partially filled pipe. The study on investigation of the 3D free surface dynamics in partially filled pipes is very limited. The method involves tinting the water white with titanium dioxide so that the water surface appears like a solid surface to image-based measuring systems. This method uses a high-resolution projector to project a stochastic pattern of light onto the water surface and uses two high-resolution cameras to capture the pattern on the deformed water surface. The 3D instantaneous water surface fluctuation can be computed from the images captured by the two cameras using a standard Digital Image Correlation algorithm. It is demonstrated that the surface dynamics parameters of turbulent flow in partially filled pipes, including surface fluctuations and surface velocity, can be measured using the projector, high-resolution cameras and DIC algorithm.


Introduction
The measurement of a liquid free surface is fundamentally important in many applications. For example, the characterization of ship wakes (Caplier et al. 2019), the measurement of ocean waves (Wanek and Wu 2006), the dynamics of sessile water drops (Chang et al. 2013), the nature of a walker's wave field (Eddi et al. 2011) and understanding of surface features of turbulent flows over a rough bed Nichols et al. 2016;. The techniques to measure the liquid surface fluctuations can be divided into two categories: intrusive methods and non-intrusive methods. Traditional intrusive methods include buoys and wave gauges (also known as wave probes). Non-contact methods are mainly based on optical systems (Moisy et al. 2009;Nichols et al. 2016), radio (Welber et al. 2016), and acoustic waves (Krynkin et al. 2014;Nichols et al. 2013;Horoshenkov et al. 2013Horoshenkov et al. , 2016. Surface measurement techniques using optical devices have been continuously developed and widely used in the last few decades (von Häfen et al. 2021;Jain et al. 2021;Nichols et al. 2020), with improvements in image resolution, acquisition speed, memory storage and decreasing transfer time. Optical measurement techniques have the advantage of being non-intrusive while allowing the comprehension of the space-time dynamics of waves over an area rather than at selected points in the case of conventional measurement systems such as wave probes.
The principle of optical techniques to measure dynamics of liquid surfaces may be broadly divided into three categories: refraction-based methods, particle illumination methods, and projection-based methods. The refraction-based techniques measure the slope of the free surface based on the intensity of the refracted light (Moisy et al. 2009;Zhang and Cox 1994;Savelsberg et al. 2006). It is a technique that is relatively easy to set up but not suitable for deep water or waves with strong curvature (Moisy et al. 2009). The particle illumination method was implemented by Turney et al. (2009), Douxchamps et al. (2005, Simonini et al. (2021) and involves seeding the surface with buoyant particles. Significant effort is required for seeding preparation, distribution of seeding and seeding illumination. The accuracy of this technique is also limited to the distribution of particles. Projection-based techniques are seedless methods conducted by projecting a pattern on the water surface. The water is usually tinted to make the water opaque so that the surface appears solid. Nichols et al. (2020) have the summarized literature on the use of colourant for surface visualization and suggested an optimum concentration of colourant. Different projection patterns have been implemented, such as a matrix of dots (Gomit et al. 2015), images of fringes (Cobelli et al. 2009), and speckle image patterns (Tsubaki et al. 2005).
Until recently, there has been no study on the water surface 3D measurement in partially filled pipes, likely due to the geometry complexity. Refraction-based techniques require a flat and transparent bed, which is not suitable for the perspex pipe used in this study. Pipe flow involves strong secondary currents and bursting motions near the free surface (Ng et al. 2021), which makes the uniform distribution of seeding particles floating on the water surface difficult. Therefore, the projection-based method is implemented in this study.
This paper is organized in the following manner: "Experimental facilities and flow conditions" section describes the experimental set-up. "DIC data processing" section includes a description of the data processing. The results of the novel DIC measurement method are presented in "Results" section and are discussed in "Discussion" section. The final conclusions are drawn in "Conclusion" section.

Pipe set-up
The work reported in this study has been carried out in a 290 mm inner diameter and 300 mm outer diameter perspex pipe (see Fig. 1) in the Hydraulics Laboratory of The University of Sheffield. The pipe is 20 m long and consists of 10 equal 2 m long sections (see Fig. 2). Individual sections were sealed with external rubber gaskets. The pipe joints were carefully machined to provide little or no disturbance to the free flow. The slope of the pipe was fixed at 1/1000 and control of the discharge from the pump was achieved with an adjustable butterfly valve (a Fisher type 8580 valve with a type 2052 actuator) upstream of the inlet tank. In Fig. 2, the labels 'a' to 'c' and 'h' to 'k' represent components of the pipe which allow control of the flow. There were several rectangular and circular entry slots cut into the pipe top to allow access for the measurement equipment. Sections 'd', 'f', and 'g' represent the measurement sections of the pipe, where the velocity field and free surface measurements take place. The measurement apparatus are described in the next section.

Measurement apparatus
The measurement of water surface fluctuation in a partially filled pipe was taken using the proposed stereoscopic DIC system and seven pairs of wave probes for validation. The DIC system used in this study is a Dantec Q400 (Dantec Dynamics) and the associated software Instra 4D. This system has been implemented in many studies for solid material deformation measurement (Mrówka et al. 2021;Gschwandl et al. 2019;Techens et al. 2020;Yan et al. 2019). As is shown in Fig. 3, the DIC system and the wave probe arrays were placed in the same pipe section. As the DIC system is a non-intrusive technique, it was placed upstream of the wave probe array to avoid any potential interference.
As is shown in Fig. 3 rightmost, a high-definition (HD) projector (NEC M403H with 1080p native resolution) was placed above the pipe projecting a 5 megapixel speckle pattern (see Fig. 4) onto the water surface and two cameras were placed just below the projector with a fixed angle in between ( ∼40 • ). The HD projector was fixed at the same position for all the tests and was focused on the water surface for each flow condition. Two high-speed CMOS cameras (Mikrotron EOSENSE 4CXP, 4 megapixel in resolution, 2336 × 1728 pixels at 563 fps frame rate, 8-bit, black-and-white) were mounted to a round bar (see Fig. 3), which allows adjusting of distance and angle in between the cameras. The round bar was attached to vertical bars, which allow the adjusting of the vertical distance of the cameras. The aperture of both cameras were always set to aperture size f8 and they were focused independently on the same patch of water surface. The sampling frequency of the DIC cameras was set to 120 Hz in order to synchronize with the HD projector which ran at 60 Hz. The overlapping field of view of the two cameras covered an area of water surface of approximately 297 mm (±31 mm) × 198 mm (±10 mm). In order to allow cameras to observe the water surface, Nichols et al. (2020) suggested at least 0.01% mass concentration of TiO2 is required for   reliable Kinect sensor measurement of surface shape. In this study, the water was tinted with titanium dioxide at a mass concentration of 0.06% for DIC measurement. The TiO2 powder was mixed in the downstream tank and was recirculated continuously by the pump. The downstream tank was agitated regularly to avoid deposition of the TiO2 powder.
A wave probe is an intrusive device that comprises two thin wires perpendicular to the water surface. Krynkin et al. (2014) analyzed the wave probe measurement error relates to the finite distance between two wires. Large error may occur when the characteristic period of surface dynamics is short. The probes were connected to standard wave probe control modules provided by Churchill Controls. The control modules energize the wave probes and measure the conductance between two wires. The instantaneous voltage between these two wires is proportional to the instantaneous water surface fluctuation. As is illustrated in Fig. 5, a wave probe calibration was completed using 14 different water depths. The water depth is proportional to the mean wave probe monitor output. This calibration procedure allows a linear regression line to be fitted empirically (dashed line in Fig. 5), which is then able to convert instantaneous voltage fluctuations to instantaneous wave fluctuations. There are 7 pairs of wave probes placed in the centre of the flume downstream of the DIC system and two wires (0.25 mm in diameter) in each pair separated 10 mm apart in the lateral direction. 10 mm distance between wires were chosen to reduce measurement error (Krynkin et al. 2014) and eliminate potential electric interference. The distances between each pair were set to be non-uniform and thus able to give 22 unique separations. The separations range from the self-paired ones (0 mm) to the furthest paired ones (500 mm). Apart from the selfpaired 0 mm separation, the shortest separation (20 mm) was placed at the downstream end of the array to reduce potential wave probes interference. To compare the water surface fluctuation simultaneously measured by two different instruments, a Dantec timing box was used. This is an interface between the control computer and the sensors, which consists of a DAQ timing and analogue data recording board (NI6211 with 16 bits resolution). The timing functions can be used to synchronize multiple sensors and ensure that the DIC images are recorded at the same time. Beside the synchronization, the box is able to record analogue signals up to 8 channels. Channels A0-A7 were used for the 7 wave probes recorded in this study.
Apart from the surface measurement apparatus, a sidelooking acoustic Doppler velocimeter (ADV) manufactured by Nortek was used to measure the 3D instantaneous velocity time series at a given spatial point under the surface. The ADV was placed 9.1 m and 14.7 m downstream the inlet tank (see Fig. 2d) to validate the uniform flow conditions. It was mounted to a modified point gauge frame and was able to adjust the depth-wise position accurate to 0.1 mm.

Flow conditions
Fourteen uniform flows were established in this study. The magnitude of the discharge was determined using an electromagnetic flow meter (Arkon MAG 910). The depth of the flow was controlled with an adjustable gate at the downstream end of the pipe to ensure uniform flow conditions throughout the measurement section. The flow depths ranged between 45.5 and 182.8 mm with bulk flow velocity from 0.30 to 0.60 m/s. The Reynolds number was between 0.84 × 10 4 and 5.05 × 10 4 as summarized in Table 1. The mean depth is the mean of the water depth measured by point gauge at the 7 streamwise locations with average spacing 1.14 m, the bulk flow velocity is calculated by U b = Q∕A and the Reynolds number is computed by where Q is the flow rate, A is the flow cross-sectional area,  is the water density, R h is the hydraulic radius, and f is the dynamic viscosity of water. Figure 6 shows the water surface captured by the upstream camera with different concentrations of TiO2. The marker 'X' in Fig. 6 corresponds to the 'X' marker in Fig. 4. In the pure water condition (Fig. 6a), no marker can be seen. As is shown in Fig. 6b and c, the marker can be barely seen and is blurry with ∼ 0.01% TiO2. The contrast and clarity is improved with higher amounts of TiO2 (see Fig. 6d, e, and f). Nichols et al. (2020) have found that at least 0.01% mass concentration of TiO2 is required to give reliable time and space averaged water depth by Kinect sensor. More than 0.01% mass concentration of TiO2 is required for this study to achieve not only time and space averaged water depth but also accurate water surface fluctuations. A compromise must be found that gives good data while negligibly affecting the surface dynamics. Flow condition 14 was running continuously and certain amount of TiO2 was gradually added into the water (from 0 to 0.04%). Figure 7 shows the power spectrum calculated from wave probe data for flow condition 14. The dashed line represents the power decay ∝ f −5 . The shape of the power spectrum and the slope of the spectrum is not affected by the increase concentration of TiO2. Titanium dioxide with mass concentration of 0.06% was used in this study and the surface dynamics were shown to not be affected by this amount of titanium dioxide. A snapshot of two camera views with 0.06% TiO2 and the synchronized wave probe voltages are illustrated in Fig. 8. The stochastic pattern can be seen clearly in both camera views while some white spots at the edges are visible which correspond to the points of direct reflection of the projector lamp.

Experimental protocol
The main experiments consist of 14 uniform flow conditions as mentioned in "Flow conditions" section. Before the turbulent flow measurements, the uniform flow condition was first examined by a side-looking acoustic Doppler velocimetry (ADV) sub-surface velocity field measurement, which will be described in detail in "Uniform flow" section. A wave probe calibration was performed before the measurement. All measurements were taken at least half an hour after the flow was established to ensure the stability of the flow (statistics of wave probe data were found to be fully developed after this time). Figure 9 shows the standard deviation water surface fluctuations measured by wave probe time series of increasing time duration as a percentage of the standard deviation of a 300 s long recording for different flows. 300 s is sufficiently long enough for the standard deviation of the surface fluctuation signals to comfortably settle to within ±1%. The surface fluctuation was measured by DIC system and wave probes synchronously for 66 s, sampled at 120 Hz and repeated 5 times for each flow. That means over 5 minutes measurement in total for each flow, which is sufficiently long enough to give convergent results. Different flow conditions have different flow depth, which means a different distance from the water surface to the camera. For each flow condition, the projection and cameras were adjusted to the focus on the water surface. A DIC system calibration was performed for each flow condition as the lens focus had been adjusted.

Calibration
Before any measurement it is necessary to perform a system calibration. A calibration target (AI-11-BMB-9× 9) from Dantec Dynamics (the black and white board shown in Fig. 10) was placed inside the pipe at the measurement area. In the Istra4D software, the appropriate type of calibration target was selected and the new calibration set-up was launched. The target was moved to each side to cover the whole measurement area and was turned in various directions and positions to identify the specific indications. Ten images of the plate in different positions and angles were used to enable each 3D calibration. The calibration plate was always fixed in the same position for the first image of the calibration to ensure the same datum and coordinate system for all flow conditions. As is shown in Fig. 10, the centre of the calibration board corresponds to (0, 0, 0) in DIC coordinates, with red and green arrows indicating streamwise and lateral directions, respectively. When cameras had registered its position, calibration parameters were saved. Uncertainties of the calibration parameters lead to measurement errors. If the value of calibration residuum was not in a good range (below 0.3-0.5) (DantecDynamics 2018), the calibration was be repeated. After some recordings of the calibration plate in various positions, the software automatically registers the nodal points of the calibration plate (see Fig. 10 green markers on the calibration board) and calculates the intrinsic (focal length, principle point and radial and

Image evaluation
After the calibration procedure, a series of images were captured by the calibrated cameras. The upper limit of measurement time was limited by the size of RAM and size of SSD/ HDD and depends on the frame rate and the image resolution. After the measurement, the evaluation was performed in the in-built Istra4D software. Images are then meshed into several facets (a small squared image region containing a set of pixels) as is shown in Fig. 11 with 25 pixels in facet size and 17 pixels in grid spacing. The correlation algorithm is based on the grey value pattern in small local neighbourhoods. The correlation algorithm can determine the shift, rotation and distortion of these facet elements (Martin and Miroslav 2021). This correlation algorithm is performed for each individual facet in the measured surface. With the known intrinsic and extrinsic parameters obtained from the calibration, the 3D displacement of each facet of the water surface can be determined.

Post processing
The raw DIC export from the Istra 4D software for each time step is the (x, y, z) coordinate of each facet relative to the initial position of the calibration plate as described in "Calibration" section. The centre of the calibration plate is 500 mm upstream of wave probe 1. Therefore the DIC xy coordinate (0, 0) means the point in the lateral centre of the flow and 500 mm upstream the wave probe 1. The raw DIC data points are non-equidistantly spaced, with some missing points in the areas with reflections or air bubbles. The first step of the DIC data pre-processing is to interpolate the Fig. 10 Illustration of the calibration in Istra4D, with the upstream camera view on the left and the downstream camera view on the right Fig. 11 Mesh the recorded images into several facets raw DIC data with the MATLAB function "scatteredInterpolant" (MATLAB 2019) into orthogonal equidistantly spaced data as is shown in Fig. 12a and b. The scatter interpolation usually does not function well in areas with missing data and leads to spikes in the interpolated surface. Tsubaki et al. (2005) observed similar protrusions in their measurement of a fluctuating free surface using a projection-based technique, which are caused by the matching error of local surface images. They suggested to be eliminated by a statistical procedure. In this study, to eliminate the spikes in the data caused by reflections and air bubbles, a 2D median filter and a 10 Hz low-pass filter were applied to the interpolated data (see Fig. 12c). Figure 13 shows the difference of the filtered data and the interpolated data over 2D space for flow 12, repeat 1 at frame 55. It shows high differences at positions like (88, −6), (100, −38), and (156, 0), where the interpolation function failed to interpolate the missing data points. Despite these areas, the filter procedure does not change the value of the z coordinate by more than 0.2 mm. The time series of z coordinate fluctuations at the xy coordinates (100, −38) before and after filtration is demonstrated in Fig. 14. The difference before and after filtering are generally below 0.2 mm, showing that spikes points can be removed robustly while not losing flow turbulence information.
The surface fluctuation is derived by subtracting the z coordinate fluctuations from the mean of the z coordinate fluctuations and then subtract the mean value. The mean surface of the z coordinate fluctuations for 14 different flows is shown in Fig. 15a. The mean of these surfaces over 2D space and the water depth measured by point gauge are illustrated in Fig. 15b. The vertical error bars presented were estimated Fig. 12 The free surface measured by DIC system for flow condition 12 repeat 1 frame 55: a raw, b interpolated, and c filtered Fig. 13 The difference between interpolated and filtered data over 2D space for flow 12 repeat 1 at frame 55 Fig. 14 The difference between interpolated and filtered data over time for flow 12 repeat 1 at spatial point (100, −38) as the maximum difference between the results obtained from the five reproducibility experiments. The maximum difference for 14 flows is 0.41 mm, which shows the repeatability of the measurement. The dashed line represents the best fit line between the DIC mean z coordinate and the point gauge measured water depth. It shows a linear relationship with gradient close to 1 between the DIC measurement and the point gauge measurement.

Results
In this section, the hydraulic conditions are first examined. After that, the comparison of DIC and wave probes measurement are made in terms of probability density, roughness height and power spectrum. Then the ability of the DIC system to recover the instantaneous free surface and surface space-time matrix is examined.

Uniform flow
As the surface fluctuations measured by DIC and wave probes were at different locations, it is important to ensure the flow is uniform in the measurement area. As is illustrated in Fig. 16a, the value of dimensionless quantities U b ∕U bf and Q∕Q f derived from experimental measurement overlap well on the theoretical graph proposed by Camp (1946), where U bf and Q f are bulk flow velocity and flow rate for pipe flowing full. Figure 16b shows the flow rates against AR h 2∕3 S 1∕2 from Manning's equation with a best fit line, where S is the slope of the pipe. The best fit line shows a best fit Manning's roughness value of 0.0099, which approximates the value of perspex (the material the pipe is made from). Furthermore, a side-looking ADV was placed 0.9 m upstream of DIC measurement area and 2.7 m downstream of the wave probes measurement area. The sub-surface velocity fluctuations were measured by it at two different streamwise positions and 15 different depthwise positions. ADV was sampled at 100 Hz and recorded by the software Vectrino Plus for a 300 s measurement duration. Details of the ADV set-up can be checked in Wu et al. (2022). The velocity profiles measured at these two positions for the highest flow, flow condition 14, are compared in Fig. 16c. Markers represent the mean streamwise velocity and error bars correspond to the standard deviation of velocity fluctuations. As is shown in Fig. 16c, the mean streamwise velocity profile and standard deviation of velocity fluctuations are comparable upstream and downstream of the measurement area. The downstream mean streamwise velocity is slightly smaller in magnitude compared to upstream one, with a small mean streamwise velocity difference within 3.95%. The difference likely arises from the intrusive wave probes in between these two ADV measurement positions. The difference is larger in the depth-wise positions 0.3 < z∕d < 0.65 and relatively smaller in the near surface and bottom positions. where is the water surface fluctuation, calculated from the instantaneous elevation series minus its mean. represents the standard deviation of surface fluctuation, also known as surface roughness height. The width of the PDF (water surface wave height) increases with the increase of flow. The width of the PDF is smaller from DIC in low flows (see Fig. 17a and b) while wider in high flows (see Fig. 17e and f).

Surface roughness
The surface roughness height (standard deviation) of surface fluctuations for all flow conditions obtained by DIC against wave probes is shown in Fig. 18 with a 1:1 linear line for reference. The surface roughness height was calculated from: (1) where T is the water surface fluctuation measurement period. The surface roughness height at a single location measured by DIC and wave probe are termed as DIC and wp , respectively. The mean DIC surface roughness height DICm is averaged over 50 streamwise locations with 2 mm spacing on average at lateral centre-line as Eq. 3 and mean wave probes surface roughness height wp m is averaged over all 7 pairs as Eq. 4: The error bars in Fig. 18 represent the maximum and minimum values obtained by DIC at different streamwise locations or different wave probes. It is evident that these two techniques show the most comparable results for medium flows. DIC shows higher surface roughness in high flows than wave probes and vice versa for low flows. Figure 19 shows the plot of the standard deviation of the surface fluctuations measured by both techniques against the ratio of water depth and pipe inner diameter d/D. Both techniques reveal a parabolic trend of the relationship between the surface roughness and the water depth.  techniques show higher spectral amplitudes for higher flows and a similar spectral pattern. It can be seen that DIC system suggests lower spectral amplitudes in low flows and higher amplitudes in high flows. Figure 21 shows three example snapshots of free surface fluctuations measured by DIC at time intervals of 0.25 s. It is evident that the surface is rougher and changes more rapidly for higher flows, which agrees with Fig. 19. A detailed instantaneous free surface fluctuation with surface gradient vectors in a 100 × 100 mm area at a time interval of 0.04 s is shown in Fig. 22. The gradient vector shows local direction and magnitude of surface slope, calculated using MATLAB "quiver" function. Figure 23 shows the space time matrices of instantaneous water surface elevation at lateral centre-line obtained by DIC system for flow conditions 1, 6, and 14. The greyscale represents the instantaneous deviation from zero. It is a visual assessment of the behaviour of the dominant features. The dashed lines correspond to the surface velocity measured by timing a floating tracer over a 3.83 m distance downstream the surface measurement area. The surface velocity is plot for every 2 s duration for better visualization with the greyscale streaks.

Discussion
The results presented in "Results" section enable a detailed exploration of the potential of DIC to resolve free surface features in partially filled pipes. The results will be discussed here.
The PDF plots in Fig. 17 compare the wave probe and DIC data and demonstrate that the DIC method tends to show higher density in small fluctuations for low flow rates, suggesting that the DIC is better able to detect small surface fluctuations than the wave probes. This can be explained by the 10 mm distance between the two wires of the wave probe. This separation distance is likely to average the surface features between two wires (Krynkin et al. 2014), and thus, they are unable to detect features smaller than the separation between wires.
As shown in Fig. 19, the standard deviation of surface fluctuations increases with the increase of flow depth, following a parabolic trend. This relationship between surface roughness and water depth was also observed by Krynkin et al. (2014), Horoshenkov et al. (2013) in rectangular open channel flows. The standard deviation measured by wave probes for the highest flow is somewhat offset from the parabolic trend. This is likely due to strong lateral motions in this flow affecting the wave probe measurement. It is also possible that there are wake effects from upstream wave probes on downstream probes for this high flow condition ( U b = 0.6 m/s). However,  has empirically tested that the intrusive wave probe wires have a negligible impact on water surface roughness for a flow with U b = 0.6 m/s in a rectangular open channel. As is shown in Fig. 24, the standard deviation of surface standard deviations calculated from different streamwise locations for DIC and wave probes is generally below 0.05, which means a small variation in wave heights between different streamwise locations. However, large variance appears between different wave probes at the highest two flow conditions, supporting the occurrence of wave probe interference in high flows. Furthermore, wave probes show larger surface roughness than the DIC when the pipe is less than half filled, while relatively smaller surface roughness than DIC when more than half filled. The difference between results from DIC and wave probes is not surprising as the DIC measures in small areas while the wave probes measure between two separated points (as discussed in "Probability density function" section). In addition, the vertical standard deviation over a 2D space measured by DIC (see Fig. 25) shows that there is some variance in both streamwise and lateral directions. It can be observed that strong lateral variation occurs in the high flow condition (Q = 28 L/s). This is likely due to secondary currents (Ng et al. 2021) at this depth ratio and further study is required.
As shown in Fig. 20, the difference in the spectral amplitudes in low and high flows can also be explained by the separation between the two wires of the wave probe as mentioned in "Probability density function" section. The slope of the power spectrum is comparable in the 5-10 Hz range and suggests that the frequency density of the energy decay is proportional to the power of −5. The −5 spectral slope, known as a power function law, was also observed by Dolcetti (2016), Tani and Fujita (2020) in rectangular open channel flows.
All three examples of instantaneous free surface fluctuation in Fig. 21 show the advection of the water surface pattern from upstream to downstream (left to right in the figure). Some advection of the surface pattern in the lateral direction can also be observed in moderate (Fig. 22b) and high flow regimes (Fig. 22c). It also shows that the size of the surface pattern is larger and advection speed is faster in higher flows. Nichols (2014) has applied quantitative U-level analysis on the surface data measured by laser-induced fluorescence (LIF) in open channel flows to estimate structure length and advection velocity. It was found that both structure length and advection velocity increase with the increase of water depth, which agrees with the observations in this study.
There are three types of motion that can be seen from the space time matrices in Fig. 23: turbulence features travelling at the speed of the mean surface flow velocity (corresponding to the red dashed lines in Fig. 23), and travelling and receding waves (corresponding to the zig-zag patterns in Fig. 23a and b). In the same spatial area, about 4 periods of oscillation can be observed for flow condition 1 (Fig. 23a), about 2 periods for flow condition 6 ( Fig. 23b) and an oscillation is hard to observe in flow condition 14 (Fig. 23c). Similar space time matrix patterns were observed by Nichols  Tani and Fujita (2020), Fujita et al. (2011), andFujita (2017). The trains of features fluctuating over space and time can be explained by the oscillation theory proposed by Nichols et al. (2016). The spatial frequency of the oscillation was shown to decrease with the increase of flow depth and increase of surface roughness, agreeing with the results here.

Conclusion
The 3D deformation of water surfaces in a partially filled pipe was measured using a DIC technique in this study. The basic principle of this technique is to project a stochastic speckle pattern onto the tinted opaque water surface and to image the surface from two cameras above. The displacement of each small facet can be evaluated and thus the fluctuations across the whole surface can be determined. This technique is able to measure surface fluctuations in a 2D area with approximate size 297 × 198 mm at high frequency (120 Hz). Larger surface area may be achieved with a different projector and camera set-up, though the measurement resolution will be limited by camera field of view and pixel density. Nichols et al. (2020) concluded that a TiO2 concentration over 0.01% is required for reliable optical surface measurement using a Kinect sensor. It was shown that a TiO2 concentration above 0.01% can affect the fluid properties (such as viscosity) and >1% is more significantly affected. In the present study, a mass concentration of 0.06% TiO2 was found adequate to make the water opaque enough for optical measurement via DIC, while not significantly changing the turbulent surface dynamics in partially filled pipes. The measured surfaces from DIC and wave probes are comparable in terms of the PDF, standard deviation and power spectrum, validating the ability of the DIC system to measure free surface dynamics at multiple locations on the surface. This technique is more suitable for surface measurement in partially filled pipes than the refraction-based technique, as it does not need to take the curvature of the pipe into account. Furthermore, it is a seedless method, which costs less and saves effort in seeding preparation. The shortcoming of this technique is that it requires opaque water, and thus the simultaneous optical measurement of sub-surface properties such as the velocity field is difficult, though other methods (such as ADV) can be used. In addition, this technique may be difficult to implement in the field due to the need for opaque water, though application may be found in highly turbid flows such as wastewater. Coupling this novel quantification of free surface dynamics with the advanced theory of Nichols (2014), Dolcetti (2016) may enable cost-effective remote assessment of flow conditions in a range of applications.