∎

11institutetext: M. Ratz, C. Cierpka 22institutetext: Technische Universität
Ilmenau, Institute of Thermodynamics and Fluid Mechanics, Germany
33institutetext: D. Fiorini, A. Simonini, M.A. Mendez (✉) 44institutetext: von
Karman Institute for Fluid Dynamics, Belgium
44email<EMAIL_ADDRESS>
# Analysis of an unsteady quasi-capillary channel flow with Time Resolved PIV
and RBF-based super resolution

Manuel Ratz  Domenico Fiorini  Alessia Simonini  Christian Cierpka and *Miguel
A. Mendez

###### Abstract

We investigate the interface dynamics in an unsteady quasi-capillary channel
flow. The configuration consists of a liquid column that moves along a
vertical 2D channel, open to the atmosphere and driven by a controlled
pressure head. Both advancing and receding contact lines were analyzed to test
the validity of classic models for dynamic wetting and to study the flow field
near the interface. The operating conditions are characterized by a large
acceleration, thus dominated by inertia. The shape of the moving meniscus was
retrieved using Laser-Induced Fluorescence (LIF)-based image processing while
the flow field near was analyzed via Time-Resolved Particle Image Velocimetry
(TR-PIV). The TR-PIV measurements were enhanced in the post-processing, using
a combination of Proper Orthogonal Decomposition (POD) and Radial Basis
Functions (RBF) to achieve super-resolution of the velocity field. Large
counter-rotating vortices were observed, and their evolution was monitored in
terms of the maximum intensity of the $Q$-field.
The results show that classic contact angle models based on interface velocity
cannot describe the evolution of the contact angle at a macroscopic scale.
Moreover, the impact of the interface dynamics on the flow field is
considerable and extends to several capillary lengths below the interface.

###### Keywords:

Unsteady quasi-capillary ChannelTime Resolved PIV RBF super resolutionDynamic
contact angles Dynamic wetting

##  1 Introduction

The dynamics of a gas-liquid interface moving along a wall play an essential
role in many industrial applications such as wetting and dewetting processes
(e.g. inkjet printing, see Eral et al. (2013)), coating applications (e.g.
slot die coating, see Kistler and Schweizer (1997)) or capillary driven flows
(e.g. liquid absorption in porous media, see Ma et al. (2012)).

(a)

(b)

Figure 1: Schematic sketch of the stream lines in proximity of the contact
line for a wetting configuration ($\theta<90^{o}$) in case of an advancing
contact line (left) and a receding contact line (right). The flow pattern on
the left is known as _stream line split injection_ , the flow pattern on the
right is known as _stream line split ejection_. Images adapted from Fuentes
and Cerro (2005).

Sufficiently close to the contact line, i.e. the intersection of the gas-
liquid interface and the solid bounding walls, capillary forces play a leading
role in the dynamics of the interface motion and thus the liquid spreading. At
a macroscopic level, capillary forces depend on the curvature of the gas-
liquid interface, but their computation requires boundary conditions at the
wall: these are usually given in terms of _contact angle_ $\Theta$, that is
the angle between the gas-liquid interface and the solid wall at the contact
line (cf. Figure 1).

Because the contact angle cannot be predicted in the framework of continuous
mechanics, one must resort to non-equilibrium statistical mechanics, _ad hoc_
assumptions in the continuous framework (E. B. Dussan V, 1976; Cox, 1998;
Voinov, 1977) or empirical laws that link the contact angle to macroscopic
quantities such as the contact line velocity (e.g. Hoffman (1975); Kistler
(1993)) or acceleration (e.g. Bian et al. (2003); Ting and Perlin (1995)).
While these approaches have been successful in ‘quasi-steady’ conditions and
very viscous fluids (Quéré, 1997; Wu et al., 2017), their generalization to
accelerating contact lines and low-viscosity fluids remain an open research
topic for both capillary or quasi-capillary flows (see Fiorini et al. (2021);
Quéré (1997); Shardt et al. (2014); Willmott et al. (2020)) and impacting
droplets (Bartolo et al. (2005); Wang and Fang (2020)).

A well known challenge in the modelling of flows near contact lines is the
singularity that arises at the contact line when combining the notion of _no-
slip_ , i.e. zero velocity at the wall, with the notion of _dynamic contact
line_ (Dussan V. and Davis, 1974; Huh and Scriven, 1971). A possible
resolution to this paradox is the splitting streamline pattern proposed by Huh
and Scriven (1971). This pattern can be modeled as a _creeping flow_ obeying a
bi-harmonic equation in the stream-function and is schematically illustrated
in Figure 1 for advancing and receding contact lines. In an advancing
configuration (Figure 1a), the fluid is pushed from the bulk of the flow
towards the wall, while in a receding configuration (Figure 1b), the opposite
pattern is expected. In both cases, the splitting streamlines pattern produces
a rolling motion near the contact line and has been observed in various
experimental campaigns (Chen et al. (1997); Fuentes and Cerro (2005); Nasarek
et al. (2008); Zimmerman et al. (2010)). In an open channel flow, such as
those encountered in capillary-rise problems (Washburn, 1921; Levine et al.,
1976), the presence of the rolling motion collides with the assumption of
fully developed flows. It is thus relevant to analyze how far from the
interface one might expect to see the impact of capillary forces near the
contact line (Savelski et al., 1995; Nasarek et al., 2008).

(a)

(b)

Figure 2: Fig (a) Sketch of the channel showing the area of the TR-PIV
measurements. The Field Of View (FOV) and Region of Interest (ROI) are
discussed in detail in Section 3.2. Fig (b) shows a photography of the
experimental set-up with the camera (1), the channel (2), the reservoir (3),
the release valve (4), the pressure ports for the injection (5 and 6) and
discharge (7).

The scope of this work is twofold. On the one hand, we aim at analyzing the
well-posedness of a correlation linking the macroscopic contact angle to the
kinematics of the contact line in inertia dominated conditions, i.e. in the
presence of large accelerations and large velocities. On the other hand, we
aim to analyze the flow field in the proximity of the contact line and measure
up to which distance from the interface the rolling motion is present. Both
aspects are analyzed considering a 2D, quasi-capillary channel (i.e. with
width $\delta\approx 2l_{c}$, with $l_{c}=\sqrt{\sigma/\rho g}$ the capillary
length, $\sigma$ the surface tension, $\rho$ the liquid density and $g$ the
gravitational acceleration). The experimental set-up and measurement
conditions are recalled in Section 2. The configuration of interest is a
narrow rectangular channel in which water moves along a solid wall forming two
dynamic contact angles. The flow is sustained by a controlled pressure head
that reproduces both advancing and receding contact lines.

The shape of the interface was tracked using Laser-Induced Fluorescence
(LIF)-based image visualization (Planquart et al. (2002)) and edge detection,
while the contact line was characterized via an inverse method by fitting a
simple model for the interface shape (Fiorini et al. (2021)). The set-up for
the LIF-based flow visualization and the related image processing is described
in Section 3.1. The velocity field was characterized via Time-Resolved
Particle Image Velocimetry (TR-PIV). Vortices near the contact line were
visualized and their strength quantified in terms of $Q$-fields.

To accurately compute the spatial derivatives required in the $Q$-field
definition, we propose a super-resolution method that combines Proper
Orthogonal Decomposition (POD) and Radial Basis Function (RBF) regression
(Karri et al. (2009); Raben et al. (2012)). The set-up for the TR-PIV
measurements and the post processing of the velocity fields is described in
Section 3.2. The results of both experimental campaigns are presented in
Section 4 while conclusions and perspectives are collected in Section 5.

##  2 Experimental Set-up and Test Cases

A sketch of the experimental set-up for the measurements is shown in Figure
2a. This consists of a rectangular channel with width
$\delta=$5\text{\,}\mathrm{mm}$$, in the figure depth
$L=$250\text{\,}\mathrm{mm}$$ and height $H=$150\text{\,}\mathrm{mm}$$.

This channel is open to atmosphere at the top and connected to a pressurized
chamber on the bottom. The chamber sustains the flow via a (time varying)
gauge pressure $\Delta p(t)$. To reproduce an advancing contact line (i.e.
interface moving upwards, along the $y$-axis), a pressure step of amplitude
$\Delta P$ is introduced by the opening of a fast electronic valve connected
to a pressure line. To reproduce a receding contact line (i.e. interface
moving downward), the same procedure is followed in reverse: the chamber is
initially pressurized, i.e. $\Delta p(0)>0$ and the release valve is opened to
the atmosphere.

In the current set-up, the opening of the valve is carried out manually.
Therefore, the experimental conditions for the receding contact line can
slightly differ between experiments. Nevertheless, the time evolution of the
pressure inside the tank is monitored with a pressure gauge and can thus be
linked to the dynamics of the interface motion as described in Section 2.1.
For the advancing contact angle problem, the fast opening of the valve (which
occurs linearly in about $0.2$ s) makes the evolution of the pressure in the
channel more repeatable, as shown in Figure 3.

A picture of the experimental setup is shown in Figure 2b. The high speed
camera is placed perpendicularly to the channel’s cross section, which is
illuminated with a laser sheet from the right side. The camera and the laser
sources changed between the LIF-visualization and the PIV experiments, as
described in Sections 3.1 and 3.2 respectively.

Figure 3: Temporal evolution of the gauge pressure in the pressure tank with a
pressure difference $\Delta$P = $1500\text{\,}\mathrm{Pa}$ for two different
experiments.

The injection/removal of air during the pressurization/depressurization used
in the advancing/receding experiments was performed on two sides of the
chamber to ensure the symmetry of the flow and minimize the entry effects.
Together with the camera (1), the picture shows the channel (2) and the
pressurized chamber (3), the fast opening valve (4), the two pressure ports (5
and 6) and the pressure gauge (7).

The walls of the channel were made in acrylic glass and the experiments were
carried out using demineralized water (density
$\rho=$997\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$$, dynamic
viscosity $\mu=$0.001\text{\,}\mathrm{Pa}\text{\,}\mathrm{s}$$, and surface
tension $\sigma=$0.072\text{\,}\mathrm{mN}\text{\,}{\mathrm{m}}^{-1}$$). The
liquid was introduced into the facility via a lateral port, and all
experiments were carried out on an initially dry surface. The influence of the
surface pre-conditioning was tested by repeating the advancing experiments
after intervals of approximately $15$ minutes. Since no noticable difference
was observed in the interface dynamics (nor in the contact evolution) between
the first and the following experiments, we conclude that this time is long
enough to have dry surfaces.

The static contact angle was measured using the same LIF-based image
processing approach described in Section 3.1 on a quiescent meniscus. This
approach is analogous to the Meniscus Profile Method (MPM) as proposed by
Petrov and Sedev (1993). The interface was positioned in the camera’s field of
view by gently rising the chamber pressure, and the measurement was taken
after waiting long enough to have a static interface. The static contact angle
was found to be $\Theta_{0}=33^{o}\pm 2^{o}$.

Finally, we emphasize that the need for splitting the experimental campaigns
into a high-speed LIF and a TR-PIV stems from the different objectives and the
different technical constraints for the two measurements. The LIF-based
visualization aims to analyze the motion of the gas-liquid interface during
its rise and oscillation, while the TR-PIV seeks to explore the flow field
during the passage of the interface. Consequently, the LIF-based visualization
required videos of longer duration, lower frame rate and larger field of view
than the TR-PIV campaign. Before proceeding with the description of the
measurement set-up, it is worth elaborating on a simple attempt to model the
response of the liquid column in the next section; this also allows to better
understand the test conditions investigated.

###  2.1 Modeling and Parameter Definition

Models of the capillary rise have been proposed by Levine et al. (1976);
Washburn (1921) and were used for the current campaign to have a preliminary
evaluation of the interface displacement as a function of the imposed pressure
head.

Denoting as $m$ the mass of the water column in motion and as $h(t)$ the
liquid height as a function of time, the force balance acting at the inlet of
the channel sets

$\displaystyle\frac{\text{d}}{\text{dt}}\left(m\dot{h}\right)=F_{g}+F_{\Delta
p}+F_{\text{s}}+F_{\text{v}}\,$ (1)

where $\dot{h}$ is the time derivative of the column height (thus the mean
velocity of the flow), $F_{g}$ is the force due to gravity, $F_{\Delta p}$ is
the force induced by the gauge pressure in the chamber, $F_{\text{s}}$ is the
force due to capillary forces at the interface, and $F_{\text{v}}$ is the
viscous force exerted at the channel walls.

Assuming that the contact line moves at the same velocity as the average
liquid height in the channel, the inertial term on the left hand side of (1)
can be further expanded as:

$\displaystyle\frac{\text{d}}{\text{dt}}\left(m\dot{h}\right)=\frac{\text{d}}{\text{dt}}\left(\rho
L\delta h\dot{h}\right)=\rho L\delta\left(\dot{h}^{2}+\ddot{h}h\right)\,,$ (2)

where the dots denote differentiation in time.

The first two terms on the right hand side do not require particular
assumption expect for the uniformity of the pressure head and the liquid
height. These are:

$\displaystyle F_{g}$ $\displaystyle=-mg=-\rho\delta Lhg$ (3) $\displaystyle
F_{\Delta p}$ $\displaystyle=\Delta p(t)A_{\text{proj}}=\Delta p(t)L\delta\,,$
(4)

where $A_{\text{proj}}$ is the channel cross-section and the liquid mass is
$m=\rho\delta Lh$.

The surface tension contribution $F_{\text{s}}$ in equation (1), depends on
the shape of the meniscus, which in turn depends on the (unknown) contact
angle. Denoting as $f(x)$ the interface shape with respect to the interface
location at the center of the channel, any curvature of the interface in the
span-wise direction is neglected and the Young-Laplace equation gives the
capillary induced pressure difference at the interface:

$\displaystyle\Delta
p_{\text{s}}=\frac{\sigma}{\delta}\int_{-\delta/2}^{\delta/2}\kappa(x)\text{d}x,\quad\kappa(x)=\frac{f^{\prime\prime}}{\left(1+{f^{\prime}}^{2}\right)^{\frac{3}{2}}}$
(5)

where $\kappa$ is the interface curvature, averaged across the channel width,
$f$ denotes the interface shape and the prime denotes spatial derivatives (see
also section 3.1. The computation of this term can be done once a model of the
interface shape is considered. The model used in this work is described in
detail in Section 3.1, and the capillary contribution to equation (1) is
$F_{\text{s}}=p_{\text{s}}A_{\text{proj}}$.

Finally, the contribution of the viscous forces $F_{\text{v}}$ in equation (1)
depends on the shape of the stream-wise velocity profile near the wall.
Assuming that this profile remains parabolic at all times, a shear stress can
be estimated:

$\displaystyle\tau$ $\displaystyle=-\mu\partial_{x}u=\frac{-4\mu
u_{\text{max}}}{\delta}=-\frac{12\mu}{\delta}\dot{h}$ (6)

where the maximum velocity is linked to the mean velocity ($\dot{h}$) using
the assumption of parabolic velocity profile. The contribution of the wall
shear stress in equation (1) is thus $F_{\text{visc}}=\tau
A_{\text{wet}}=-h\dot{h}{12L\mu}/{\delta}$, with $A_{\text{wet}}$ being the
wetted surface of the channel walls.

Given a model for the interface shape $f(x)$ (which depends on the evolution
of the dynamic contact angle $\Theta(t)$), and given the time-varying pressure
in the chamber $\Delta p(t)$, it is possible to integrate equation (1) and
compute the evolution of the liquid height under the aforementioned
assumption. While the presented model is oversimplified (it does not account,
among other things, for the pressure losses at the channel entrance), it is
possible to analyze the relative impact of the various terms to the interface
dynamics. Moreover, this model was used to link, at least approximately, the
pressure in the chamber to the velocity and acceleration of the contact line
during the various testing conditions.

###  2.2 Experimental Conditions and Test Cases

In the LIF-based flow visualization the experiments were conducted in four
conditions, namely three advancing cases and one descending case. The
advancing test cases were carried out setting $\Delta p(0)=0$ and $\Delta
P$=1000, 1200 and $1500\text{\,}\mathrm{Pa}$. These led to $\Delta
p(t\rightarrow\infty)$=660, 825 and $985\text{\,}\mathrm{Pa}$ respectively.
The descending case was carried out with $\Delta
p(0)=$1200\text{\,}\mathrm{Pa}$$ and $\Delta p(t\rightarrow\infty)=0$. The
relevant experimental parameters are listed in Table 1. In these experiments,
the camera position was not varied: this was placed at the top of the channel,
with the upper edge of the image at approximately $15\text{\,}\mathrm{mm}$
from the upper edge.

All experiments were repeated three times. Since no appreciable differences
were observed in the interface dynamics, we report in this article only the
result of a single run per test case. Concerning the descending test case, the
choice of presenting only one case is due to the limited impact observed by
the initial step $\Delta p(0)$ on the evolution of the interface shape: during
the descent, the interface tends to reach a constant velocity after a certain
distance and this velocity is fairly independent of the initial pressure step.
The chosen value is the one that allows for observing the interface over a
large field of view, providing enough space for the video acquisition at
multiple heights.

Table 1: Experiments in the LIF-based visualization. LIF | $\Delta p(0)$ | $\Delta P$ | $\Delta p(t\rightarrow\infty)$ | camera pos.
---|---|---|---|---
Test Case 1 | 0 | 1000 | 660 | 0
Test Case 2 | 0 | 1200 | 825 | 0
Test Case 3 | 0 | 1500 | 985 | 0
Test Case 4 | 1200 | - | 0 | 0
Table 2: Experiments in the TR-PIV campaign. TR-PIV | $\Delta p(0)$ | $\Delta P$ | $\Delta p(t\rightarrow\infty)$ | camera pos
---|---|---|---|---
Test Case 1 | 0 | 1500 | 985 | 0
Test Case 2 | 200 | 1200 | 985 | 0
Test Case 3 | 1200 | - | 0 | 0
Test Case 4 | 1200 | - | 0 | 40
Test Case 5 | 1200 | - | 0 | 120

For the PIV measurements, two rising cases and three descending cases were
considered. The experimental parameters are listed in Table 2. The ‘Test Case
1’ in the PIV campaign corresponds to the conditions of ’Test Case 3’ in the
LIF-visualization campaign. In the Test Case 2, the acceleration was reduced
by reducing the pressure step setting $\Delta p(0)=200$ Pa. The other test
cases are in the same conditions as ‘Test Case 4’ in the LIF-visualization
campaign, and the only difference is in the position of the camera which is
moved downward by $40\text{\,}\mathrm{mm}$ and $120\text{\,}\mathrm{mm}$ with
respect to the reference value. As for the LIF visualization experiments, also
the TR-PIV experiments were repeated three times and no significant
variability was observed.

##  3 Measurement Techniques

###  3.1 LIF Visualization: Set-up and Processing

For the LIF measurements, Rhodamine B from Sigma Aldrich is dissolved in the
water. Images were captured with a Phantom v2012 high-speed camera which has a
resolution of 1200$\;\times\;$800 px. The particles were illuminated with a
continuous Stabilite 2017 ion laser system with a wavelength of
$515\text{\,}\mathrm{nm}$. The power of the laser was set to
$1.9\text{\,}\mathrm{W}$. An objective lens with a focal length of
$105\text{\,}\mathrm{mm}$ was used. The camera was placed at a larger distance
from the channel compared to the PIV experiments, to allow the tracking of the
interface over a large vertical distance. The achieved optical magnification
was 30 px/mm. Images were acquired in single frame mode, with a frame rate of
$500\text{\,}\mathrm{Hz}$ and an exposure time of
$20\text{\,}\mathrm{\SIUnitSymbolMicro s}$. In total, 1500 images were
acquired, resulting in a measurement duration of $3\text{\,}\mathrm{s}$.

(a)

(b)

Figure 4: Fig (a) Raw image of the meniscus obtained from the experiments. Fig
(b) Meniscus after non-local means denoising and high-pass filtering.

A zoom view of the LIF-based visualization of the dynamic meniscus is shown in
Figure 4. The images were processed using the non-local means denoising by
Buades et al. (2011) and then high-pass filtered to highlight the gas-liquid
interface. The result of this step, for the image in Figure 4a, is shown in
Figure 4b. Only the right side of the symmetric meniscus is shown. The
interface is detected by calculating the vertical gradient of the image
intensity, similarly to Mendez et al. (2016, 2018), by searching for the peak
in the intensity gradient in each image column. At the end of this process,
outliers were removed using a local median filter.

A simple model for the interface shape is then fitted to the detected points.
This model was heuristically adapted from the analytical solution of a
meniscus in stationary conditions and reads:

$f(x)=\text{cosh}\left(\frac{x^{a}}{b}\right)-1,$ (7)

where $y=f(x)$ is the vertical position of the interface with respect to the
interface at the center of the image (see Figure 5) and $a,b\in\mathbb{R}$ are
constants to be fit at each time step.

Figure 5: Location of the reference frame moving with the meniscus, and with
respect to which the model in equation (7) is defined.

The fitting is performed using the minimize function from Python’s Scipy
library (Virtanen et al., 2020) to minimize the $l_{2}$ norm of the prediction
error and the time evolution of each parameter is further smoothed using a
low-pass filter to eliminate outliers.

Equation 7 is used to (1) compute the contact angle by solving
$\text{ctg}(\Theta)=f^{\prime}(\delta/2)$, with ′ denoting the spatial
derivative, (2) to compute the location of the contact line as
$h(\delta/2,t)=f(\delta/2)+h(t)$, with $h(t)$ the liquid height at the center
of the channel and (3) compute the contact line velocity as
$u_{c}(t)=\dot{h}(\delta/2,t)$.

###  3.2 TR-PIV: Set-up

For the PIV measurements, the fluid was seeded with Red Fluorescent Polymer
Microspheres from Thermo Fisher Scientific with a diameter of
$12\text{\,}\mathrm{\SIUnitSymbolMicro m}$. The flow was observed with a TR-
PIV system from Dantec Dynamics. The particles were illuminated with the laser
(DM40-527-DH Nd:YLF from Photonics Industries) with a wavelength of
$527\text{\,}\mathrm{nm}$ and a maximum pulse energy of
$40\text{\,}\mathrm{mJ}$ at a frequency of $1\text{\,}\mathrm{kHz}$. Images
were recorded with a SpeedSense Ethernet M310 camera having a resolution of
1280$\times$800 px (at a maximum frame rate of $3260\text{\,}\mathrm{Hz}$).
The acquisition frequency of both laser laser and camera, operating in _single
frame mode_ , was and $1.2\text{\,}\mathrm{kHz}$. An objective lens with a
focal length of $105\text{\,}\mathrm{mm}$ was used to get an optical
magnification of 60 px/mm.

Figure 6 shows two snapshots from the PIV acquisition _before_ (a,c) and
_after_ (b,d) the image pre-processing, which was carried out using the POD-
based background removal by Mendez et al. (2017). Figure 6a shows a snapshot
in which the interface is not within the field of view; in this case the
illumination and the seeding was fairly homogeneous. Figure 6c shows a
snapshot in which the interface is in the field of view; in this case the
interface shape and the light refraction hinder the interface detection and
the velocimetry in its proximity.

Figure 6: PIV snapshots before and after the pre-processing. The interface is
not in the field of view in figure (a) while it is in the field of view in
(c). Figures (b) and (d) result from the eigenbackground removal from (a) and
(c) respectively.

This is why the PIV interrogation was carried out in a cropped region of the
FOV moving with the interface and having its upper boundary (edge 4-3 in
Figure 2a) approximately $1.5\text{\,}\mathrm{mm}$ below the expected liquid
level. This distance was measured in several of the snapshots for which the
interface is visible and was kept constant by a tracking algorithm. In
particular, the Region of Interest (ROI) in which the PIV was carried out, is
displaced for each image pair $k$, at time $t_{k}$, using the mean velocity
computed along the bottom edge 1-2 (denoted as $y_{1,2}$), i.e.
$\overline{v}_{k}=\int_{-\delta/2}^{\delta/2}v(x,y_{1,2},t_{k})dx$. The
vertical displacement of the FOV is thus $\Delta x=\overline{v}\Delta t$, with
$\Delta t$ the sampling frequency of the PIV. The velocity of the ROI was
found to closely match the mean velocity of the interface $h(t)$.

In the ROI, the PIV interrogation was carried out with the OpenPIV software by
Liberzon et al. (2020), using windows of large aspect ratio, justified by the
fact that velocity component $v$ (along $y$) is much larger than the velocity
component $u$ (along $x$). The rectangular interrogation windows allows for a
better sampling of the flow in the direction where the largest gradients are
expected, without deteriorating the signal-to-noise ratio in the cross-
correlation maps.

Adaptive iterative multigrid interrogation (Scarano and Riethmuller, 2000) was
used to calculate the velocity fields. For the advancing contact line
experiments, the initial window size is set to 256$\times 64$ px and reduced
over two steps to 64$\times 16$ px. For the receding contact line experiments,
the initial size is 64$\times$64 px, and reduced to 24$\times$24 px in two
iterations. Outliers were removed based on the signal-to-noise ratio measured
in terms of peak to standard deviation in the cross-correlation map.

###  3.3 Superresolution of PIV fields

The post-processing of the PIV fields aimed at increasing the resolution of
the PIV data and enable accurate computation of $Q-$fields. The approach
consist in using a Radial Basis Function (RBFs, see Fornberg and Flyer (2015))
expansion of the velocity field to compute derivatives analytically. A similar
approach has been used for robust interpolation (Casa and Krueger, 2013), to
compute derivatives near walls (Karri et al., 2009) and compute pressure
fields from PIV (Sperotto et al., 2021). In this work, we use it to compute
the $Q-$fields analytically and thus accurately measure the intensity of the
vortices near the gas-liquid interface.

To reduce the computational cost of the RBF regression, we combine it with a
classic Proper Orthogonal Decomposition (POD). The main idea behind the POD-
RBF-super-resolution strategy is to perform RBF regression of the spatial and
temporal structures of the decomposition, then use these to reconstruct high-
resolution POD modes and finally rebuild high-resolution velocity fields. The
procedure is herein briefly summarized. Let

$\vec{u}(\textbf{x}_{0},t_{0})=\sum_{i=1}^{R}\sigma_{i}\vec{\phi_{i}}(\textbf{x}_{0})\varphi_{i}(t_{0})\,,$
(8)

denote the truncated POD of the velocity field on an initial spatial grid
$\mathbf{x}_{0}$ containing $n_{s}$ points and a temporal grid $t_{0}$ with
$n_{t}$ points. Here $R\ll n_{t}$ denotes the truncation index, $\sigma_{i}$
is the amplitude of the i-th POD mode with spatial structure
$\vec{\phi_{i}}(\textbf{x}_{0})$ and temporal structure $\varphi_{i}(t_{0})$.

The details of the POD computation can be found elsewhere (e.g. Mendez et al.
(2019, 2020a)) and are omitted here. We denote as

$\vec{\phi}_{i}(\textbf{x}_{0})=\sum_{j=1}^{n_{\phi}}w^{\phi}_{ij}\,\gamma^{\phi}_{j}(\mathbf{x}_{0}|\mathbf{x}_{i,j},\Sigma_{\phi_{j}}),,$
(9)

and

${\varphi}_{i}(t_{0})=\sum_{j=1}^{n_{\varphi}}w^{\varphi}_{i,j}\,\gamma^{\varphi}_{j}(t_{0}|t_{j},\Sigma_{\varphi_{j}})\,$
(10)

the RBF expansion of the spatial and temporal structures respectively.
Therefore, $w^{\phi}_{i,j}$ and $w^{\varphi}_{i,j}$ are the set of weights
defining the regression functions,
$\gamma^{\phi}_{j}(\mathbf{x}_{0}|\mathbf{x}_{j},\Sigma_{\phi_{j}})$ are the
$n_{\phi}$ radial basis functions in space and
$\gamma^{\varphi}_{j}(t_{0}|t_{j},\Sigma_{\varphi_{j}})$ are the $n_{\varphi}$
radial basis functions in time. These have collocation points $\mathbf{x}_{j}$
and $t_{j}$ and shape parameters $\Sigma_{\phi_{j}}$ and
$\Sigma_{\varphi_{j}}$. Note that the coefficients $w^{\phi}_{i,j}$ must be
defined as vectors, i.e. $w^{\phi}_{i,j}=(w^{\phi}_{u,i,j},w^{\phi}_{v,i,j})$
since the same radial basis $\gamma^{\phi}_{j}$ is used for both components
$\vec{\phi}=(\phi_{u},\phi_{v})$. However, because every snapshot is reshaped
as a column vector, regardless of whether this collects a vector or a scalar
field, we keep the same notation as for the interpolation of the temporal
structures. We consider Gaussian RBFs both in space and in time. The spatial
RBFs $\gamma^{\phi}_{j}(\mathbf{x}_{0}|\mathbf{x}_{j},\Sigma_{\phi_{j}})$ are
defined as

$\gamma^{\phi}_{j}=\exp\left(-(\mathbf{x}_{0}-\mathbf{x}_{j})^{T}\Sigma^{-1}_{\phi_{j}}(\mathbf{x}_{0}-\mathbf{x}_{j})\right)$
(11)

with
$\Sigma^{-1}_{\phi_{j}}=\text{diag}(1/2\sigma^{2}_{x},1/2\sigma^{2}_{y})$,
while the temporal RBFs
$\gamma^{\varphi}_{j}(t_{0}|t_{j},\Sigma_{\varphi_{j}})$ are defined as

$\gamma^{\varphi}_{j}=\mbox{exp}\left(-(t_{0}-t_{j})^{2}/(2\sigma^{2}_{t})\right)\,$
(12)

with $\Sigma_{\varphi_{j}}=1/2\sigma^{2}_{t}$.

For both space and time regressions, the collocation point and the shape
parameters are defined a priori and the regression is solved once the weights
are identified. Reshaping all bases functions as columns of matrices
$\Gamma_{\phi}(\mathbf{x}_{0})\in\mathbb{R}^{n_{s}\times n_{\phi}}$ and
$\Gamma_{\varphi}(t_{0})\in\mathbb{R}^{n_{t}\times n_{\varphi}}$, collecting
all the unknown weights into column vectors
$\mathbf{w}^{\phi}\in\mathbb{R}^{n_{\phi}}$,
$\mathbf{w}^{\varphi}\in\mathbb{R}^{n_{\phi}}$, and collecting all the entries
of $\vec{\phi}_{i}(\textbf{x}_{0})$ and ${\varphi}_{i}(t_{0})$ into column
vectors $\mathbf{\phi}_{0}\in\mathbb{R}^{n_{s}}$ and
$\mathbf{\varphi}_{0}\in\mathbb{R}^{n_{t}}$, the weights are the solution of
classic least square problems. Using a Tikhonov regularization to mitigate the
risks of over-fitting (Mendez et al., 2020b), the solution is:

$\mathbf{A}\mathbf{x}=\mathbf{b}\quad\rightarrow\mathbf{x}=\Bigl{(}\mathbf{A}^{T}\mathbf{A}+\alpha\mathbf{I}\Bigr{)}^{-1}\mathbf{A}^{T}\mathbf{b},$
(13)

where $\mathbf{A}=\Gamma_{\phi}(\mathbf{x}_{0})$,
$\mathbf{x}=\mathbf{w}^{\phi}$ and $\mathbf{b}=\mathbf{\phi}_{i}$ for the
regression in space and $\mathbf{A}=\Gamma_{\varphi}(t_{0})$,
$\mathbf{x}=\mathbf{w}^{\varphi}$ and $\mathbf{b}=\mathbf{\varphi}_{i}$ for
the regression in time. The regularization parameter $\alpha$ is taken as
$\epsilon\,||\mathbf{A}||_{F}$, with $\epsilon=1e-11$ and $||\,||_{F}$ the
Frobenious norm. Note that the inverse in (13) can be pre-computed and the
system solved for all the modes in one single step.

Given the weights, equations (9) and (10) can be used on _any_ arbitrary mesh
grid $\mathbf{x}$ and _any_ time discretization $t$. Then, the resulting high
resolution mode can be introduced in the expansion (8).

The analytic description of the flow allows for the computation of derivatives
by simply replacing the RBF with the required (analytic) derivatives. For
example, $\partial_{x}u$ can be computed by replacing
$\vec{\phi}_{i}(\textbf{x})$ with $\partial_{x}\vec{\phi}_{i}(\textbf{x})$ in
(8) and this can be computed by replacing
$\gamma^{\phi}_{j}(\mathbf{x}|\mathbf{x}_{j},\Sigma_{\phi_{j}})$ with
$\partial_{x}\gamma^{\phi}_{j}(\mathbf{x}|\mathbf{x}_{j},\Sigma_{\phi_{j}})$
in (9). These derivatives are analytically available because
$\partial_{x}\gamma^{\phi}_{j}$ is analytically available by differentiating
(11).

Finally, using the available derivatives it is possible to compute the
strength of vortices in the PIV fields in terms of $Q$-field, which for a 2D
incompressible flow reads (Hunt et al., 1988)

$\begin{split}Q&=\frac{1}{2}(||\Omega||_{F}^{2}-||S||_{F}^{2})\\\
&=-\frac{1}{2}\left(\left({\partial_{x}u}\right)^{2}+2{\partial_{y}u}\,\partial_{x}v+\left({\partial_{y}v}\right)^{2}\right)\end{split}\,$
(14)

where $S=1/2(\nabla V+\nabla V^{T})$ and $\Omega=1/2(\nabla V-\nabla V^{T})$
the symmetric and anti symmetric portions of the velocity gradient tensor.

(a)

(b)

(c)

Figure 7: LIF results for Test Case 1. Fig (a) corresponds to the interface
above the velocity field in Figure 10d. The three images are marked by the
three red squares in Figure 9e. The titles recall the $Ca$ and $G$ numbers for
each snapshot as well as the time $t$ at which the meniscus came into the FOV
of the camera.

(a)

(b)

(c)

Figure 8: LIF results for Test Case 4. In each snapshot, the values of $Ca$
and $G$ are given as well as the time $t$ at which the meniscus came into the
FOV of the camera.

##  4 Results

###  4.1 LIF Visualization Results

We begin by illustrating the results of the interface detection via high speed
LIF-based visualization. The results of the repeated experiments yielded
similar results, thus the results always show the case of the first
repetition. Figure 9 shows three snapshots of Test Case 1, each taken at
different points of the oscillation. Figure 8 shows three snapshots of Test
Case 4. In each figure, the text in the image recalls the dimensionless
conditions characterizing each snapshot, namely the dimensionless velocity of
the contact line in terms of capillary number $Ca=\mu u_{c}/\sigma$ and
dimensionless acceleration $G=a/g$ as well as the time passed after the
meniscus entered the FOV. The signs of $Ca$ and $G$ follow the coordinate
system defined in Fig. 2, i.e. a positive sign indicates an upwards motion or
an acceleration upwards. The results of the interface fitting, following
equation (7) is shown using a dashed line.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 9: LIF results of the oscillating interface for different test cases.
Fig (a), (c) and (e) temporal evolution of the mean interface height inside
the channel for Test Case 1, 2 and 3 respectively. The dashed line represents
the height taken from the images, the dotted line the predicted height from
the model in Section 2.1. Fig (b), (d) and (f) show the measured contact angle
as a function of capillary number $Ca=u_{c}\mu/\sigma$. The markers are
colored according to the acceleration number $G$. The capillary number at
$t=$0\text{\,}\mathrm{s}$$ is $Ca=5.3\times 10^{-3}$, but the axis is limited
in the range $[-1.25,1.25]\times 10^{-3}$ to better visualize the final region
at $Ca\approx 0$. The red squares in Fig. (f) correspond to images in Figure 7
(c.f the times label in the captions of that figure).

Despite the widely different conditions in terms of velocity and acceleration,
the shape of the interface is well described by the model in equation (7). The
interface remains symmetric with respect to the vertical axis, while the
contact angle in dynamic conditions largely varies between the rising and the
descending conditions.

For Test Case 1 to 3, Figure 10 plots on the left the evolution of the liquid
height as a function of time (blue dashed lines) together with the prediction
of the simple model in section 2.1. These test cases differ in the amplitude
of the pressure step introduced, thus on the overall level of acceleration.
Despite the simplifications involved, the model correctly predicts the
interface’s position and velocity in the first $0.5$ s, i.e. at the end of the
first cycle. Afterwards, as the interface oscillates around its equilibrium
conditions, the model prediction loses accuracy.

Among the various terms in equation (1), it was found that the inertial
contribution (left-hand side term) was the dominant one, followed by the
viscous and the surface tension terms. Given the good performances of the
simple model in equation (7) in representing the shape of the interface and
thus the accurate prediction of the capillary contribution, it is clear that
the main limitation of the model arises from the modelling of the viscous
term. Although this was not within the scope of this work, this result should
be considered for future works on unsteady capillary or quasi capillary
channels (e.g. see Washburn (1921); Levine et al. (1976)).

Finally, on the right side of Figure 10, maps for every test case the
evolution of the contact angle $\Theta$ as a function of the capillary number
$Ca$ is shown. The markers are coloured by the level of acceleration $G$ in
each point. The range of contact angles spanned within one experiment is
$\Theta\in[20,110^{o}]$. The largest values are produced during the first rise
of the interface when the acceleration is also the largest ($G=-0.24$). The
range of capillary number observed during the experiments is
$Ca\in[-1.25,5.30]\times 10^{-3}$ and is approximately the same in all test
cases. Large capillary numbers are not shown in the plots to focus on the
regions with small capillary numbers. During the descending phase, the
capillary number varies between $Ca=-1.25\times 10^{-3}$ to $Ca\approx 0$; yet
the contact angle remains approximately constant and equal to $\Theta\approx
22^{o}$. On the contrary, during the final phase of the experiment, when the
interface oscillation are vanishing and the contact line is almost at rest,
the contact angle exhibits the most prominent variation (in the range
$\Theta\in[20,80]^{o}$. In this phase, also the acceleration is approximately
constant. While these results suggest that the history of the contact line
dynamic might play a role in predicting the contact angle, it also brings the
question of whether it is possible to use classic contact line models of the
form $\Theta=f(Ca,G)$ for low viscous fluids and in the presence of
acceleration.

###  4.2 TR-PIV Experiment Results

(a)

(b)

(c)

(d)

(e)

Figure 10: PIV results for the oscillating interface. Fig (a) - (b)
dimensionless acceleration ($G$) and maximum $Q$-field as a function of the
capillary number $Ca$ for two test cases. The blue circles correspond to Test
Case 1 and the orange crosses to Test Case 2. Fig (c) - (e) $Q$-field contour
and high-pass filtered velocity field. The aspect ratio of the plots is one.

(a)

(b)

(c)

(d)

(e)

Figure 11: Results for the descending interface. The figure is the same as
Figure 10, but considering Test Case 3, 4 and 5. The labels refer to the
vertical position of the camera, where $0\text{\,}\mathrm{mm}$ is the
configuration in which the bottom boundary of the FOV is at the level of the
reservoir’s upper edge (see Figure 2).

The PIV velocity fields were processed retaining a total of $R=25$ POD modes
in equation (8). The temporal structures of the modes were approximated with
$n_{\varphi}=500$ RBFs with $\sigma_{t}=0.01$ while $n_{\phi}=900$ RBFs were
used for the spatial structures (see equation (11)). For the regression in
space, the RBFs were collocated over a regular grid of $n_{Bx}\times
n_{By}=60\times 15$ and anisotropic kernels with $\sigma_{x}=10$ and
$\sigma_{y}=15$ were selected. This allows to efficiently account for the high
aspect ratio of the flow. No re-sampling was performed in the time domain
since the available time resolution was sufficient for the scope of this work.

Figure 10 shows the results of Test Case 1 and 2, i.e. the advancing contact
line. Figure 10a and 10b show the gravitational acceleration $G$ and the
maximum of the $Q-$field representing the vortex strength within the ROI as a
function of the capillary number $Ca$. The labels in the legend refer to the
operating conditions described in Section 2.1. The capillary number in these
plots is different from the one used in the LIF campaign, i.e. it  is not
based on the velocity of the contact line but on the mean velocity computed by
zero-padding the velocity field at the walls and taking the mean of the
average velocity of the bottom five rows. This operation acts as smoothing and
serves for plotting purposes. We distinguish these quantities with an over-bar
in the labels of figure axes.

A clear trend is visible: the $Q$-field increases with decreasing mean flow
velocity.

Three instantaneous velocity fields are shown in Figures 10c \- 10e,
corresponding to the three points marked in Figure 10a and 10b. For plotting
purposes, the velocity fields are slightly high-pass filtered by using a
Gaussian filter with $\sigma_{x}=\sigma_{y}=15$ and a truncation after $2.5$
standard deviations. The quiver plot only shows every third vector along $x$
and every second one along $y$ for better visibility. The $Q$-field is also
shown in each snapshot for quantitative analysis, and the axis aspect ratio is
set to one. We recall that the upper horizontal boundary of the ROI is at
approximately $1.5\text{\,}\mathrm{mm}$ from the interface at the centre of
the channel. In each case, two counter-rotating vortices are visible close to
the walls (the one on the right rotates clockwise and the one on the left
counterclockwise). The vortices extend considerably along the vertical
direction but only over a small distance in the cross-stream direction.
Therefore, they have no remarkable influence on the velocity field at a
distance of $1\text{\,}\mathrm{mm}$ from the wall.

The initial pressure was high enough to produce a strong acceleration, as
shown in Figure 10a. While the rolling motion of the flow is still present,
the large acceleration and high velocities push the stream-line split
injection pattern and its vortices towards the wall. It is worth noticing that
the velocity fields look similar even though Figure 10c and 10d have largely
different accelerations. Referring to the expected theoretical flow field from
Figure 1a, it appears that the stream-line split injection forms a smaller
angle with the walls than what pictured in the schematic of Figure 1a. As a
result, the rolling motion expected by the viscous-capillary balance is
confined to a narrow region of the flow.

Figure 11 shows the same results for Test Case 3 to 5. The capillary and
acceleration numbers in Figure 11a and 11b have a different sign because the
interface moves downward. A similar trend as for the rising interface can be
observed, with the magnitude of the $Q$-field increasing as the velocity
decreases. However, the reader should note that a decreasing velocity, in this
case, means that the absolute value of the velocity is increasing. Another
observation is the increased strength of the $Q$-field, with the maximum value
being almost an order of magnitude above the maximum value for the rising
interface.

The vortices in Figure 11c \- 11e show a different behaviour than the ones for
the rising interface. Besides the reversed rotation, as expected from the
theoretical flow topology in Figure 1, their centres are much closer to the
centre of the channel, and the rolling motion extends much further in both the
$x$ and $y$ directions. In the case of the receding contact line, the stream-
line split ejection (see Figure 1b) forms a larger angle with the wall than
for the advancing contact line in the investigated configuration
($\Theta_{0}=33\pm 2^{o}$).

The large difference cannot be justified by the different velocity, which in
the slowest rising test cases is comparable to the fastest falling test cases,
nor by its acceleration, which is also similar in the cases in 11c and the one
in Figure 10e. Instead, the main difference between the rising and descending
test cases is in the shape of the interface (cf. Figure 8 a) for a snapshot
during the rise and Figure 8 c) for a snapshot during the fall). This impacts
the velocity flow field in the entire ROI, thus more than
$5\text{\,}\mathrm{mm}$ below the interface.

##  5 Conclusions and Perspectives

The dynamics of accelerating menisci was investigated using high-speed LIF-
based visualization and Time-Resolved PIV. The first was used to track the
evolution of the interface shape and the dynamic contact angle, while the
second was used to measure the flow field in the fluid adjacent to the
interface. Both advancing and receding contact lines were investigated for a
wetting configuration (static contact angle $\Theta_{0}=33\pm 2^{o}$)
consisting of water flowing over acrylic glass.

A simple model was used to describe the interface shape over the full range of
operating conditions, with capillary numbers in the range
$Ca\in[-1.25,5.30]\times 10^{-3}$ and contact lines acceleration in the range
$G\in[-0.24,0.07]$. This model was also used to measure the dynamic contact
angle which reaches values of $\Theta=110^{o}$ during a rise (at large
accelerations) and $\Theta=20^{o}$ during a descent. Within the entire range
of investigated conditions, it appeared impossible to relate the dynamic
contact angle to the contact line velocity and/or acceleration: during some
phases of the experiment, the contact angle remains constant over a wide range
of $Ca$ and $G$, while in others it varies significantly over small intervals
of $Ca$ and $G$. The quest for identifying a functional dependency
$\theta=f(Ca,G)$ in inertia-dominated conditions might thus not be well posed.

Concerning the velocity field underneath the interface, two counter-rotating
vortices were observed below the menisci in both advancing and receding cases.
These were detected, and their intensity quantified in terms of $Q-$field,
using a super-resolution approach combining RBFs and POD. The topology of
these vortices complies with the well known splitting streamline pattern
postulated to solve the apparent incompatibility of contact line motion versus
no-slip condition at the wall. The strength of these vortices depends on the
mean velocity and the acceleration of the flow, and significantly differs
between advancing and receding conditions. In the receding conditions, the
splitting streamline forms a much larger angle with the wall and has a much
larger impact on the flow field. Significant differences were observed between
advancing/receding wetting configurations even when the modulus of velocity
and acceleration of the contact lines were comparable. This highlights the
relevance of capillary forces near the interface, since the main difference
between the two cases was found in the interface’s shape: this appeared nearly
flat in the advancing cases and was characterized by large menisci during the
receding test cases.

## References

  * Bartolo et al. (2005) Bartolo D, Josserand C, Bonn D (2005) Retraction dynamics of aqueous drops upon impact on non-wetting surfaces. Journal of Fluid Mechanics 545(-1):329, DOI 10.1017/s0022112005007184
  * Bian et al. (2003) Bian X, Schultz WW, Perlin M (2003) Liquid slug motion and contact lines in an oscillatory capillary tube. Journal of Fluid Dynamics (under revision), http://wwwpersonal umich edu/schultz/Manuscripts/BianSchultzPerlin pdf 146:151
  * Buades et al. (2011) Buades A, Coll B, Morel J (2011) Non-local means denoising. Image Processing On Line 1, DOI 10.5201/ipol.2011.bcm˙nlm
  * Casa and Krueger (2013) Casa LDC, Krueger PS (2013) Radial basis function interpolation of unstructured, three-dimensional, volumetric particle tracking velocimetry data. Measurement Science and Technology 24(6):065304, DOI 10.1088/0957-0233/24/6/065304
  * Chen et al. (1997) Chen Q, Ramé E, Garoff S (1997) The velocity field near moving contact lines. Journal of Fluid Mechanics 337:49–66, DOI 10.1017/S0022112096004806
  * Cox (1998) Cox RG (1998) Inertial and viscous effects on dynamic contact angles. Journal of Fluid Mechanics 357:249–278, DOI 10.1017/s0022112097008112
  * Dussan V. and Davis (1974) Dussan V EB, Davis SH (1974) On the motion of a fluid-fluid interface along a solid surface. Journal of Fluid Mechanics 65(1):71–95
  * E. B. Dussan V (1976) E B Dussan V (1976) The moving contact line: the slip boundary condition. Journal of Fluid Mechanics 77(4):665–684, DOI 10.1017/s0022112076002838
  * Eral et al. (2013) Eral H, t’ Mannetje D, Oh JM (2013) Contact angle hysteresis: a review of fundamentals and applications. Colloid and polymer science 291(2):247–260
  * Fiorini et al. (2021) Fiorini D, Mendez MA, Simonini A, Seveno D, Johan S (2021) Transient dynamics of contact angles. to appear in Journal of colloid and interface science
  * Fornberg and Flyer (2015) Fornberg B, Flyer N (2015) Solving pdes with radial basis functions. Acta Numerica 24(3):215–258
  * Fuentes and Cerro (2005) Fuentes J, Cerro R (2005) Flow patterns and interfacial velocities near a moving contact line. Experiments in Fluids 38(4):503–510
  * Hoffman (1975) Hoffman RL (1975) A study of the advancing interface. i. interface shape in liquid—gas systems. Journal of Colloid and Interface Science 50(2):228 – 241, DOI 10.1016/0021-9797(75)90225-8
  * Huh and Scriven (1971) Huh C, Scriven L (1971) Hydrodynamic model of steady movement of a solid/liquid/fluid contact line. Journal of Colloid and Interface Science 35(1):85–101
  * Hunt et al. (1988) Hunt J, Wray A, Moin P (1988) Eddies, streams, and convergence zones in turbulent flows. In: Studying Turbulence Using Numerical Simulation Databases, 2. Proceedings of the 1988 Summer Program, Stanford, California, USA, June 27 - July 22
  * Karri et al. (2009) Karri S, Charonko J, Vlachos PP (2009) Robust wall gradient estimation using radial basis functions and proper orthogonal decomposition (POD) for particle image velocimetry (PIV) measured fields. Measurement Science and Technology 20(4):045401, DOI 10.1088/0957-0233/20/4/045401
  * Kistler (1993) Kistler SF (1993) Hydrodynamics of Wetting. In: Wettability, CRC Press, chap 6, pp 311 – 430
  * Kistler and Schweizer (1997) Kistler SF, Schweizer PM (1997) Liquid film coating: scientific principles and their technological implications. Springer
  * Levine et al. (1976) Levine S, Reed P, Watson E, Neale G (1976) A theory of the rate of rise of a liquid in a capillary. In: Colloid and interface science, Academic Press, pp 403–419
  * Liberzon et al. (2020) Liberzon A, Lasagna D, Aubert M, Bachant P, Käufer T, Jakirkham, Bauer A, Vodenicharski B, Dallas C, Borg J, tomerast, ranleu (2020) Openpiv/openpiv-python: added synthetic image generator
  * Ma et al. (2012) Ma K, Liontas R, Conn CA, Hirasaki GJ, Biswal SL (2012) Visualization of improved sweep with foam in heterogeneous porous media using microfluidics. Soft Matter 8(41):10669–10675
  * Mendez et al. (2016) Mendez MA, Németh L, Buchlin JM (2016) Measurement of liquid film thickness via light absorption and laser tomography. EPJ Web of Conferences 114:02072, DOI 10.1051/epjconf/201611402072
  * Mendez et al. (2017) Mendez MA, Raiola M, Masullo A, Discetti S, Ianiro A, Theunissen R, Buchlin J (2017) Pod-based background removal for particle image velocimetry. Experimental Thermal and Fluid Science 80(108):723–729
  * Mendez et al. (2018) Mendez MA, Scelzo M, Buchlin JM (2018) Multiscale modal analysis of an oscillating impinging gas jet. Experimental Thermal and Fluid Science 91:256–276, DOI 10.1016/j.expthermflusci.2017.10.032
  * Mendez et al. (2019) Mendez MA, Balabane M, Buchlin JM (2019) Multi-scale proper orthogonal decomposition of complex fluid flows. Journal of Fluid Mechanics 870:988–1036, DOI 10.1017/jfm.2019.212
  * Mendez et al. (2020a) Mendez MA, Hess D, Watz BB, Buchlin JM (2020a) Multiscale proper orthogonal decomposition (mPOD) of TR-PIV data—a case study on stationary and transient cylinder wake flows. Measurement Science and Technology 31(9):094014, DOI 10.1088/1361-6501/ab82be
  * Mendez et al. (2020b) Mendez MA, Pino F, Fiore M (2020b) Machine learning for fluid mechanics: Challenges, opportunities and perspectives. In: Optimization Methods for Computational Fluid Dynamics, VKI Lecture Series, von Karman Institute
  * Nasarek et al. (2008) Nasarek R, Wereley S, Stephan P (2008) Flow Field Measurements Near a Moving Meniscus of a Capillary Flow With Micro Particle Image Velocimetry (µPIV). In: Proceedings of the Sixth International ASME Conference on Nanochannels, Microchannels and Minichannels, Darmstadt, Germany, June 23-25
  * Petrov and Sedev (1993) Petrov JG, Sedev RV (1993) Quasi-static shape of the steady moving meniscus—comparison of the apparent and extrapolated dynamic contact angles. Colloids and Surfaces A: Physicochemical and Engineering Aspects 74(2-3):233–241
  * Planquart et al. (2002) Planquart P, Riethmuller M, Joly A, Detry J (2002) Real time optical detection and characterization of water model free surfaces - application to the study of the mold. In: 4th European Continuous Casting Conference
  * Quéré (1997) Quéré D (1997) Inertial capillarity. EPL (Europhysics Letters) 39(5):533
  * Raben et al. (2012) Raben SG, Charonko JJ, Vlachos PP (2012) Adaptive gappy proper orthogonal decomposition for particle image velocimetry data reconstruction. Measurement Science and Technology 23(2):025303
  * Savelski et al. (1995) Savelski M, Shetty S, Kolb W, Cerro R (1995) Flow patterns associated with the steady movement of a solid/liquid/fluid contact line. Journal of Colloid and Interface Science 176(1):117–127
  * Scarano and Riethmuller (2000) Scarano F, Riethmuller M (2000) Pod-based background removal for particle image velocimetry. Experiments in Fluids 29(1):51–60
  * Shardt et al. (2014) Shardt O, Waghmare PR, Derksen J, Mitra SK (2014) Inertial rise in short capillary tubes. Rsc Advances 4(28):14781–14785
  * Sperotto et al. (2021) Sperotto P, Pieraccini S, Mendez MA (2021) A meshless method to compute pressure fields from image velocimetry. to appear in Measurement Science and Technology 2112.12752
  * Ting and Perlin (1995) Ting CL, Perlin M (1995) Boundary conditions in the vicinity of the contact line at a vertically oscillating upright plate: an experimental investigation. Journal of Fluid Mechanics 295(-1):263, DOI 10.1017/s0022112095001960
  * Virtanen et al. (2020) Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, Burovski E, Peterson P, Weckesser W, Bright J, van der Walt SJ, Brett M, Wilson J, Millman KJ, Mayorov N, Nelson ARJ, Jones E, Kern R, Larson E, Carey CJ, Polat I, Feng Y, Moore EW, VanderPlas J, Laxalde D, Perktold J, Cimrman R, Henriksen I, Quintero EA, Harris CR, Archibald AM, Ribeiro AH, Pedregosa F, van Mulbregt P, SciPy 10 Contributors (2020) SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods 17:261–272, DOI 10.1038/s41592-019-0686-2
  * Voinov (1977) Voinov OV (1977) Hydrodynamics of wetting. Fluid Dynamics 11(5):714–721, DOI 10.1007/bf01012963
  * Wang and Fang (2020) Wang F, Fang T (2020) Retraction dynamics of water droplets after impacting upon solid surfaces from hydrophilic to superhydrophobic. Physical Review Fluids 5(3):033604, DOI 10.1103/physrevfluids.5.033604
  * Washburn (1921) Washburn EW (1921) The dynamics of capillary flow. Physical Review 17(5):273–283
  * Willmott et al. (2020) Willmott GR, Briole A, Szczepaniak F (2020) Inertial capillary uptake of drops. Physical Review E 101(4):043109
  * Wu et al. (2017) Wu P, Nikolov AD, Wasan DT (2017) Capillary rise: Validity of the dynamic contact angle models. Langmuir 33(32):7862–7872, DOI 10.1021/acs.langmuir.7b01762
  * Zimmerman et al. (2010) Zimmerman J, Weislogel M, Tretheway D (2010) Micro-PIV Measurements near a Moving Contact Line. In: Proceedings of the ASME 2010 International Mechanical Engineering Congress & Exposition, Vancouver, British Columbia, Canada, November 12-18


in $\mathbb{D}_{M}((V_{\ell}\otimes\underline{U}\otimes
V_{r}^{*})^{\mathfrak{m}})$. The result now follows since the left-hand side
vanishes when acting on
$C^{\infty}(A_{\textup{reg}};(V_{\ell}\otimes\underline{U}\otimes
V_{r}^{*})^{M})$ by Corollary 6.4.

Finally, replacing the role of
$\widetilde{\mathcal{D}}_{i}^{\sigma_{\ell},\underline{\tau},\sigma_{r}}$ by
$\mathcal{D}_{i}^{\sigma_{\ell},\underline{\tau},\sigma_{r}}$ gives the result
when $\widehat{\kappa}$ in $\mathcal{A}$ is replaced by $\kappa$. ∎

For $n=2$ we obtain the following special case of Theorem 6.12.

######  Corollary 6.13.

The classical dynamical reflection term $\mathcal{A}$ (see (6.10)) acts as
zero on $(V_{\ell}\otimes U_{1}\otimes U_{2}\otimes V_{r}^{*})^{M}$,

(6.15) $\mathcal{A}|_{(V_{\ell}\otimes U_{1}\otimes U_{2}\otimes
V_{r}^{*})^{M}}=0.$

The equation (6.15) are called the classical dynamical reflection equation for
$\widehat{\kappa}$ relative to the pair $(r^{+},r^{-})$.

For $n=3$ we have

######  Corollary 6.14.

The classical dynamical reflection term $\mathcal{A}$ (see (6.10)) and the
mixed classical dynamical Yang-Baxter terms
$\textup{MYB}(s)\in\mathcal{R}\otimes(U(\mathfrak{g})^{\otimes 3})^{M}$ (see
(6.11)) satisfy for $1\leq i<j\leq 3$,

(6.16) $\bigl{(}\mathcal{A}_{ij}+\textup{MYB}(k)\bigr{)}|_{(V_{\ell}\otimes
U_{1}\otimes U_{2}\otimes U_{3}\otimes V_{r}^{*})^{M}}=0$

where $\\{k\\}=\\{1,2,3\\}\setminus\\{i,j\\}$ and $\textup{MYB}(k)$ acts in
the natural way on the tensor components $U_{1}\otimes U_{2}\otimes U_{3}$
within $(V_{\ell}\otimes U_{1}\otimes U_{2}\otimes U_{3}\otimes
V_{r}^{*})^{M}$.

The three equations (6.16) are called the coupled classical dynamical Yang-
Baxter-reflection equations for the triple $(r^{+},r^{-},\widehat{\kappa})$.

The classical dynamical reflection type equation (6.15) does not guarantee
that $\mathcal{A}_{ij}$ in (6.16) vanishes when acting on $(V_{\ell}\otimes
U_{1}\otimes U_{2}\otimes U_{3}\otimes V_{r}^{*})^{M}$. Hence in general the
reflection term in the coupled classical dynamical Yang-Baxter-reflection
equations (6.16) cannot be decoupled from the Yang-Baxter type terms.

######  Definition 6.15.

We call the equations (6.12) for $n>3$ the higher order coupled classical
dynamical Yang-Baxter-reflection equations for the triple
$(r^{+},r^{-},\widehat{\kappa})$.

######  Remark 6.16.

We expect that more generally,

(6.17)
$\mathcal{A}_{ij}+\sum_{k=1}^{i-1}\textup{MYB}(1)_{kij}+\sum_{k=i+1}^{j-1}\textup{MYB}(2)_{ikj}+\sum_{k=j+1}^{n}\textup{MYB}(3)_{ijk}$

vanishes when acting on $(V_{\ell}\otimes\underline{U}\otimes
V_{r}^{*})^{\mathfrak{m}}$ for all $n\geq 2$. This is the case when $G$ is
real split, see [20, Thm. 6.26] and Remark 4.19. In the real split case we
have $\mathfrak{m}=0$ and hence we conclude that in this case (6.17) vanishes
identically on $V_{\ell}\otimes\underline{U}\otimes V_{r}^{*}$ for all $n\geq
2$. These equations decouple, and are equivalent to the equations

(6.18) $\begin{split}\mathcal{A}_{12}&=0\qquad\hbox{ in
}\,\,\mathcal{R}\otimes\textup{End}(V_{\ell}\otimes U_{1}\otimes U_{2}\otimes
V_{r}^{*}),\\\ \textup{MYB}(k)&=0\qquad\hbox{ in
}\,\,\mathcal{R}\otimes\textup{End}(V_{\ell}\otimes U_{1}\otimes U_{2}\otimes
U_{3}\otimes V_{r}^{*})\,\,\,\hbox{ for }\,\,\,k=1,2,3,\end{split}$

which are the classical dynamical reflection equation for $\widehat{\kappa}$
relative to $(r^{+},r^{-})$ and the mixed classical dynamical Yang-Baxter
equation for $(r^{+},r^{-})$ from [20, Thm. 6.31].

##  7\. Example: $G=\textup{SU}(p,r)$.

Write $e_{ij}^{(m)}\in\textup{GL}(m;\mathbb{C})$ for the matrix unit with a
one at entry $(i,j)$ and zeros everywhere else. Write
$I_{m}\in\textup{GL}(m;\mathbb{C})$ for the unit matrix, and
$I_{m}^{-}:=\sum_{i=1}^{m}e_{i,m+1-i}^{(m)}\in\textup{GL}(m;\mathbb{C})$ for
the antidiagonal $m\times m$-matrix with ones on the antidiagonal. Fix $1\leq
r\leq p$ and write

$C:=\left(\begin{matrix}\frac{1}{\sqrt{2}}I_{r}&0&-\frac{1}{\sqrt{2}}I_{r}^{-}\\\
0&I_{p-r}&0\\\
\frac{1}{\sqrt{2}}I_{r}^{-}&0&\frac{1}{\sqrt{2}}I_{r}\end{matrix}\right)\in\textup{GL}(p+r;\mathbb{C})$

with $0$ the zero matrix of the appropriate size, and consider
$G:=\\{C^{-1}gC\,\,|\,\,g\in\textup{SU}(p,r)\\}$. Then

$\mathfrak{g}_{\mathbb{R}}=\\{X\in\mathfrak{sl}(p+r;\mathbb{C})\,\,|\,\,X^{\ast}J_{p,r}+J_{p,r}X=0\\}$

with $X^{*}=\overline{X}^{T}$ the adjoint of $X$ and

$J_{p,r}:=\left(\begin{matrix}0&0&-I_{r}^{-}\\\ 0&I_{p-r}&0\\\
-I_{r}^{-}&0&0\end{matrix}\right).$

The formula $\theta_{\mathbb{R}}(X):=J_{p,r}XJ_{p,r}$
($X\in\mathfrak{g}_{\mathbb{R}}$) defines a Cartan involution of
$\mathfrak{g}_{\mathbb{R}}$, and the diagonal matrices in
$\mathfrak{g}_{\mathbb{R}}$ form a $\theta_{\mathbb{R}}$-stable maximally
nonabelian Cartan subalgebra $\mathfrak{h}_{\mathbb{R}}$ of
$\mathfrak{g}_{\mathbb{R}}$ with associated
$\mathfrak{p}_{\mathbb{R}}$-component

$\mathfrak{a}_{\mathbb{R}}=\bigoplus_{j=1}^{r}\mathbb{R}(e_{jj}^{(p+r)}-e_{j^{\prime}j^{\prime}}^{(p+r)})\quad\hbox{
with }\quad j^{\prime}:=p+r+1-j.$

The root system $R\subset\mathfrak{h}^{*}$ of
$\mathfrak{g}=\mathfrak{sl}(p+r;\mathbb{C})$ is
$R=\\{\epsilon_{i}-\epsilon_{j}\,\,|\,\,1\leq i\not=j\leq p+r\\}$ with
$\epsilon_{i}(h):=\lambda_{i}$ for
$h=\sum_{j=1}^{p+r}\lambda_{j}e_{jj}^{(p+r)}\in\mathfrak{h}$. Take
$R^{+}=\\{\epsilon_{i}-\epsilon_{j}\,\,|\,\,1\leq i<j\leq p+r\\}$ as the set
of positive roots. We will write use the shorthand notation
$\alpha_{ij}:=\epsilon_{i}-\epsilon_{j}$ for the roots. Write
$\alpha_{i}:=\alpha_{i,i+1}$ ($1\leq i\leq p+q-1$) for the associated basis
elements of $R$, then the involution ${}^{t}\theta\in\textup{Aut}(R)$ is
determined by

$\theta(\alpha_{i})=\begin{cases}-\alpha_{i^{\prime}}\qquad&(1\leq i\leq
r),\\\ \alpha_{i}\qquad&(r+1\leq i<p),\\\ -\alpha_{i^{\prime}}\qquad&(p\leq
i<p+r-1),\end{cases}$

which in turn determines the Satake diagram of $G$.

Write $\mathfrak{a}_{\mathbb{R}}^{*}=\bigoplus_{j=1}^{r}\mathbb{R}f_{j}$ with
$f_{j}\in\mathfrak{a}_{\mathbb{R}}^{*}$ defined by
$f_{j}(e_{ii}^{(p+r)}-e_{i^{\prime}i^{\prime}}^{(p+r)}):=\delta_{i,j}$. The
restricted root system $\Sigma$ is a root system of type $\textup{BC}_{r}$
with the positive restricted roots given by

$\Sigma^{+}=\\{f_{i}\pm f_{j}\\}_{1\leq i<j\leq
r}\cup\\{f_{j},2f_{j}\\}_{1\leq j\leq r}.$

In fact, for $1\leq i<j\leq r$ and $1\leq k\leq r$ we have

$\begin{split}R_{0}&=\\{\alpha_{\ell m}\\}_{r+1\leq\ell\not=m\leq p},\\\
R_{f_{i}-f_{j}}&=\\{\alpha_{ij},\alpha_{j^{\prime}i^{\prime}}\\},\qquad\quad
R_{f_{i}+f_{j}}=\\{\alpha_{ij^{\prime}},\alpha_{ji^{\prime}}\\},\\\
R_{f_{k}}&=\\{\alpha_{k\ell},\alpha_{\ell k^{\prime}}\\}_{r+1\leq\ell\leq
p},\qquad R_{2f_{k}}=\\{\alpha_{kk^{\prime}}\\},\end{split}$

(for $r=1$ there are no restricted roots of the form $f_{i}\pm f_{j}$
($i\not=j$), and hence the middle line drops out). Write
$\xi_{i}:=\xi_{f_{i}}$ for the multiplicative character of $A$ associated to
$f_{i}\in\mathfrak{a}_{\mathbb{R}}^{*}$. The dependence on $A$ of the local
factors $r^{\pm},\kappa$ of the asymptotic boundary KZB operators
$\widetilde{\mathcal{D}}_{i;M}^{\sigma_{\ell},\underline{\tau},\sigma_{r}}$
($1\leq i\leq n$), as well as of the potential
$V^{\sigma_{\ell},\underline{\tau},\sigma_{r}}$ of the Schrödinger operator
$\widetilde{H}_{M}^{\sigma_{\ell},\underline{\tau},\sigma_{r}}$, are then
described in terms of the multiplicative characters $\xi_{i}\xi_{j}^{\pm 1}$
($1\leq i<j\leq r$) and $\xi_{k}$ ($1\leq k\leq r$), which are the characters
associated to the indivisible positive restricted roots of $\Sigma$. The
equivariance of the operators is with respect to the hyperoctahedral group
$W\simeq S_{r}\ltimes(\pm 1)^{r}$.

Recall that $\mathfrak{k}=\mathfrak{m}\oplus\mathfrak{q}$ with
$\mathfrak{q}=\bigoplus_{\alpha\in R^{+}\setminus
R_{0}^{+}}\mathbb{C}y_{\alpha}$. To further describe the local factors of
$\widetilde{\mathcal{D}}_{i}^{\sigma_{\ell},\underline{\tau},\sigma_{r}}$ and
the potential $V^{\sigma_{\ell},\underline{\tau},\sigma_{r}}$ of
$\widetilde{H}^{\sigma_{\ell},\underline{\tau},\sigma_{r}}$ one needs the
explicit description of $\mathfrak{m}$ and of the standard basis elements
$y_{\alpha}$ ($\alpha\in R^{+}\setminus R_{0}^{+}$) of $\mathfrak{q}$. We have

$\mathfrak{m}=\mathfrak{t}\oplus\bigoplus_{r+1\leq i\not=j\leq
p}\mathbb{C}e_{ij}^{(p+r)}$

with $\mathfrak{t}$ the $(+1)$-eigenspace of $\theta|_{\mathfrak{h}}$, which
is explicitly described by

$\mathfrak{t}=\\{\sum_{i=1}^{r}\lambda_{i}(e_{ii}^{(p+r)}+e_{i^{\prime}i^{\prime}}^{(p+r)})+\sum_{j=r+1}^{p}\mu_{j}e_{jj}^{(p+r)}\,\,|\,\,2\sum_{i=1}^{r}\lambda_{i}+\sum_{j=r+1}^{p}\mu_{j}=0\\}.$

The elements $y_{\alpha}$ ($\alpha\in R^{+}\setminus R_{0}^{+}$) are
explicitly given by

$\begin{split}y_{\alpha_{ij}}&=e_{ij}^{(p+r)}+e_{i^{\prime}j^{\prime}}^{(p+r)},\qquad
y_{\alpha_{j^{\prime}i^{\prime}}}=e_{ji}^{(p+r)}+e_{j^{\prime}i^{\prime}}^{(p+r)},\\\
y_{\alpha_{ij^{\prime}}}&=e_{ij^{\prime}}^{(p+r)}+e_{i^{\prime}j}^{(p+r)},\qquad
y_{\alpha_{ji^{\prime}}}=e_{ji^{\prime}}^{(p+r)}+e_{j^{\prime}i}^{(p+r)},\\\
y_{\alpha_{k\ell}}&=e_{k\ell}^{(p+r)}-e_{k^{\prime}\ell}^{(p+r)},\qquad
y_{\alpha_{\ell k^{\prime}}}=e_{\ell k^{\prime}}^{(p+r)}-e_{\ell
k}^{(p+r)},\\\
y_{\alpha_{kk^{\prime}}}&=e_{kk^{\prime}}^{(p+r)}+e_{k^{\prime}k}^{(p+r)}\end{split}$

for $1\leq i<j\leq r$, $1\leq k\leq r$ and $r+1\leq\ell\leq p$ (the
$1^{\textup{st}}$-$4^{\textup{th}}$ line describes $y_{\alpha}$ for the roots
$\alpha\in R_{f_{i}-f_{j}}^{+}$, $\alpha\in R_{f_{i}+f_{j}}^{+}$, $\alpha\in
R_{f_{k}}^{+}$ and $\alpha\in R_{2f_{k}}^{+}$ respectively).

Finally, since
$\mathfrak{s}(\mathfrak{gl}(p)\times\mathfrak{gl}(r))=\textup{Ad}(C)\mathfrak{k}$
the basis elements of $\mathfrak{k}$ can be expressed in terms of the standard
basis elements of $\mathfrak{s}(\mathfrak{gl}(p)\times\mathfrak{gl}(r))$, and
vice versa. For instance, the center of $\mathfrak{k}$ is spanned by

$\textup{Ad}(C^{-1})\left(\begin{matrix}\frac{1}{p}I_{p}&0\\\
0&-\frac{1}{r}I_{r}\end{matrix}\right)=\left(\begin{matrix}\bigl{(}\frac{1}{2p}-\frac{1}{2r}\bigr{)}I_{r}&0&-\bigl{(}\frac{1}{2p}+\frac{1}{2r}\bigr{)}I_{r}^{-}\\\
0&\frac{1}{p}I_{p-r}&0\\\
-\bigl{(}\frac{1}{2p}+\frac{1}{2r}\bigr{)}I_{r}^{-}&0&\bigl{(}\frac{1}{2p}-\frac{1}{2r}\bigr{)}I_{r}\end{matrix}\right),$

which can be expressed as

$\left(\bigl{(}\frac{1}{2p}-\frac{1}{2r}\bigr{)}\sum_{i=1}^{r}(e_{ii}^{(p+r)}+e_{i^{\prime}i^{\prime}}^{(p+r)})+\frac{1}{p}\sum_{j=r+1}^{p}e_{jj}^{(p+r)}\right)-\bigl{(}\frac{1}{2p}+\frac{1}{2r}\bigr{)}\sum_{i=1}^{r}y_{\alpha_{ii^{\prime}}},$

where the first term lies in $\mathfrak{t}$ and the second term lies in
$\mathfrak{q}$.

## References

  * [1] D. Bernard, On the Wess-Zumino-Witten model on the torus, Nucl. Phys. B303, 77–93 (1988).
  * [2] J. Cardy, Boundary Conformal Field Theory, Encyclopedia of Math. Phys. (Elsevier, 2006), arXiv:hep-th/0411189
  * [3] W. Casselman, D. Miličić, Asymptotic behavior of matrix coefficients of admissible representations, Duke Math. J. 94 (1982), 869–930.
  * [4] J. Dixmier, Enveloping algebras, Graduate Studies in Math. 11, Amer. Math. Soc.
  * [5] P. Etingof, F. Latour, The dynamical Yang-Baxter equation, representation theory, and quantum integrable systems, Oxford Lecture Series in Math. and its Appl. 29, Oxford University Press 2005.
  * [6] P. Etingof, O. Schiffmann, Twisted traces of intertwiners for Kac-Moody algebras and classical dynamical $r$-matrices corresponding to Belavin-Drinfeld triples, Math. Res. Lett. 6 (1999), 593–612.
  * [7] P. Etingof, O. Schiffmann, A. Varchenko, Traces of intertwiners for quantum groups and difference equations, Lett. Math. Phys. 62 (2002), 143–158
  * [8] P. Etingof, A. Varchenko, Traces of intertwiners for quantum groups and difference equations. I, Duke Math. J. 104 (2000), 391–432.
  * [9] G. Felder, Conformal field theory and integrable systems associated to elliptic curves, Proceedings of the International Congress of Mathematicians, Vol. 1,2 (Zürich, 1994), pp. 1247–1255, Birkhauser, Basel, 1995.
  * [10] G. Felder, C. Weiczerkowski, Conformal blocks on elliptic curves and the Knizhnik-Zamolodchikov-Bernard equations, Comm. Math. Phys. 176 (1996), 133–161.
  * [11] Harish-Chandra, Spherical functions on a semisimple Lie group I, Amer. J. Math. 80 (1958), 241–310.
  * [12] G.J. Heckman, E.M. Opdam, Root systems and hypergeometric functions. I, Compositio Math. 64 (1987), 329–352.
  * [13] G.J. Heckman, H. Schlichtkrull, Harmonic analysis and special functions on symmetric spaces, Perspectives in Mathematics, Vol. 16 (1994).
  * [14] J.E. Humphreys, Introduction to Lie Algebras and Representation Theory, Graduate Texts in Math. 9, Springer Verlag, New York, 1972.
  * [15] A.W. Knapp, Lie groups beyond an introduction, Progress in Math. Vol. 140 Birkhäuser 1996.
  * [16] M.A. Olshanetsky, A.M. Perelomov, Quantum systems related to root systems, and radial parts of Laplace operators, Funct. Anal. Appl. 12 (1978), 123–128.
  * [17] N. Reshetikhin, Degenerate integrability of quantum spin Calogero-Moser systems, Lett. Math. Phys. 107 (2017), 187–200.
  * [18] N. Reshetikhin, J.V. Stokman, Classical spin Calogero-Moser chains of Harish-Chandra type, in preparation.
  * [19] O. Schiffmann, On classification of dynamical $r$-matrices, Math. Res. Lett. 5, 13–30 (1998).
  * [20] J.V. Stokman, N. Reshetikhin, $N$-point spherical functions and asymptotic boundary KZB equations, arXiv:2002.02251.
  * [21] J.V. Stokman, Mixed classical dynamical Yang-Baxter equations and the associated reflection equation, in preparation.
  * [22] G. Warner, Harmonic Analysis on Semi-Simple Lie Groups I, Springer Verlag, 1972.
  * [23] G. Warner, Harmonic Analysis on Semi-Simple Lie Groups II, Springer Verlag, 1972.
  * [24] J.-B. Zuber, CFT, BCFT, ADE and all that, in: ”Quantum symmetries in theoretical physics and mathematics” (Bariloche, 2000), 233–266, Contemp. Math., 294, Amer. Math. Soc., Providence, RI, 2002.