We consider model problems for the tear film over multiple blink cycles with heat transfer from the posterior side of the tear film. A nonlinear partial differential equation governs the film thickness on a moving domain in one space dimension and time. One end of the tear film moves in order to mimic blinking in the eye. The film thickness is coupled with the diffusion of heat from the posterior of the film, where the underlying cornea and aqueous humor are modeled as a rectangular domain. The domain of the tear film is located on one edge of the rectangle. The resulting problem is solved using the method of lines with a Chebyshev spectral method in space. Evaporation is included in the model, with end fluxes specified to compensate for the evaporation from the film. The numerical results reveal a similarity to quantitative in vivo observations of the film dynamics and measured ocular surface temperature. Periodicity in the film and temperature dynamics is explored with different flux conditions and end motions, and a transition between periodic and non-periodic solutions is analyzed.

The ocular tear film is critical for good vision and eye health. The properly functioning tear film protects the ocular surface with moisture, helps transport waste away from the ocular surface, and provides a smooth optical surface for visual function.1 This multilayer film2,3 has an anterior layer of lipids that is in contact with air4–6 and is tens of nanometers thick on average. (Here we use anterior to mean the forward direction toward the air, posterior to mean the rearward direction into the cornea, inferior to mean the direction toward the cheek, and superior is the direction toward the forehead.) The lipid layer floats upon a primarily aqueous layer1 which is in turn on a mucin-rich region at the corneal surface.7–9 The precorneal tear film (PCTF) is the tear film located directly on the cornea, while the pre-lens tear film is found on the anterior surface of a contact lens. For this paper, we focus exclusively on the PCTF.

The PCTF is a total of a few microns thick in the center of the cornea after a blink10–12 and has a considerably thicker meniscus (about 0.1 mm or more) around the lid margins13–15 where the tear film climbs the wettable part of the eyelids. This structure must be re-formed rapidly after each blink to enable vision with minimal interruption. In this paper, we study the dynamics of the tear film and the ocular surface temperature (OST) of a normal eye over several blink cycles. The term “blink cycle” is used to mean the combined periods of a single blink, in which the superior lid moves down toward the inferior lid and then returns to its original position, together with the interblink period separating two blinks. The model we develop here is an extension of the model of that employed by Li and Braun16 in which the eye was assumed to remain open during the interblink period, and the introduction given here relies heavily on that work.

Theoretical models of tear film dynamics have recently been reviewed;17 a brief review is paraphrased here from Braun et al.18 The tear film is commonly assumed to be Newtonian, and the underlying substrate (the cornea) is assumed to be flat.18,23 Mathematical studies have incorporated a variety of important effects: surface tension;19–22 polar lipid surface concentration gradients causing the Marangoni effect;23–26 evaporation;6,27,28 wettability of the corneal surface via van der Waals terms;27,29,30 motion of the eyelids in one space dimension;24–26,31–35 lid motion and a separate lipid layer;36 heat transfer posterior to the tear film;16 and the shape of the eye opening.37,38 These effects may all contribute in different regions of the eye and at different times in the blink cycle.

In his theory of the lacrimal system, Doane39 posed that significant drainage along the lid margins begins with the lids about halfway open, and that it ends up to 3 s after the lids have fully opened. New tear fluid is supplied from the lacrimal gland, which secretes the aqueous part of the tear film, and this new fluid is observed to enter the exposed tear film from the superior lid near the outer (temporal) canthus.15,39–41 The aqueous part of tears exits via the puncta, which are small holes found near the nasal canthus. Tear film models have incorporated lacrimal gland supply and punctal drainage as well, in both one-dimensional32,33 and two-dimensional37,38 models.

We briefly comment on the tear film fluid properties. A sample of tears combining all of its components is mildly shear thinning42 but weakly elastic.43 Recent measurements of meibomian lipids (meibum) alone from various mammals show that meibum may have significant elasticity44 at room temperature but that the elasticity is very small above 32 °C.45 The contribution of the lipid layer to the overall rheology and properties of tears at in vivo temperatures is yet to be understood. The apparent surface tension of the tear-air interface and other interfacial properties are affected by surface active polar lipids46,47 that are thought to occupy the lipid-aqueous interface and to be insoluble in the aqueous layer. The polar lipids can cause upward motion of the tear fluid after a blink.23,48 However, the apparent surface tension of tears appears to be affected by proteins and other tear film components.49,50 In this paper, we simplify the action of the lipid layer to that of a limiting case of strong surface activity from an insoluble surfactant that results in the so-called uniform stretching limit or USL.24,31,32

Models of tear film that fix the temperature of the substrate predict that the temperature of the tear film surface increases slightly after a blink;27,28 however, ocular surface temperature measurements show that cooling occurs, typically about 1 or 2 °C (discussed below). In this paper, we focus on the fluid dynamics of the tear film during the blink cycle and its interaction with heat transfer inside the eye posterior to the tear film. Li and Braun16 examined three different models for the heat distribution posterior to the tear film during the interblink only; the treatments of heat transfer used there were inspired by models of Ajaev and Homsy.51–54 The substrate is considered either thick or thin depending on the value of

$y^{\prime }_c$
yc⁠, where
$y^{\prime }_c$
yc
specifies the size of the heat conduction domain. For
$y^{\prime }_c$
yc
comparable to the tear film thickness, the substrate was considered thin. For
$y^{\prime }_c$
yc
on the order of a few millimeters, the substrate was considered thick. Li and Braun16 showed that only with the thick substrate case can the tear film model explain the observed cooling rates and final temperatures for the OST observed experimentally. In this paper, we extend this model to compute dynamics during the blink and the interblink, that is, for the whole blink cycle.

Measurement of the OST has been reviewed recently.55 Non-contact temperature measurement of OST using infrared thermometry was pioneered by Mapstone;56–58 this approach appears to be the most successful method for determining the temperature of the anterior surface of the tear film. This method has been used to measure anterior surface temperatures for both the PCTF59,60 and the pre-lens tear film.61 

In this paper, we compare our results with the results of Efron et al.59 because of the quantitative reporting and interpretation of their results. They found that the location of the average minimum temperature of the cornea, based on 21 subjects, was slightly inferior to (below) the geometric center of the cornea (GCC). They speculated that this was because the eyelids heat the eye and the GCC is superior to the geometric center of the palpebral fissure. The temperature of the GCC 2s after a blink was found to be 34.3 ± 0.7 °C on average. The horizontal temperature variation across the ocular surface was fit well with a parabola. The temperature is higher at the periphery of the cornea, and toward the lid margins; the temperature of the limbus was 0.45 °C higher than at the GCC. An earlier study62 found that this temperature difference was 0.6 °C.

Efron et al.59 also measured the rate of cooling of the ocular surface; they found a rate of 0.033 ± 0.024 °C/s for the first 15s after a blink. They also found that if a subject had a slower rate of cooling and that if the OST equilibrated, then a subject was more likely to be able to refrain from blinking for an extended period. The data from their experiments suggested that the corneal temperature could vary by 0.78 °C on average for each subject, and we will find that including this temperature variation is necessary for understanding ocular surface temperature variation for thin film modeling.

Scott63,64 developed finite element (FE) models of heat transfer in the eye that were intended to help decide whether radiative heating led to cataracts in subjects employed as glass blowers. Scott developed63 a detailed axisymmetric FE model of heat transfer within the eye and to the external environment beyond the anterior of the eye. Scott64 then incorporated radiative heating from the extreme environment of furnaces used in glass-blowing to examine the effects within the eye. Other finite element models of heat transfer in the eye65,66 found steady-state OST values that agree with experimentally determined ranges. Steady-state finite element models with conductive and convective heat transfer in the aqueous humor have been developed in two67,68 and three69 dimensions with similar conclusions; the location of the coolest OST and the overall heat transfer through the anterior chamber is not significantly affected by the distorted temperature field. Approximate analytical models have been developed as well.70,71 Siggers and Ethier72 reviewed flow in the aqueous humor concerning the transport of heat and proteins and the hydration of the cornea. A major point of the paper by Li and Braun16 is that the necessary modeling needed to approximate the surface dynamics of the tear film is substantially less complicated, and a simple rectangular geometry posterior to the tear film was sufficient to capture OST dynamics. Their model suggested an initial cooling rate of 0.058 °C/s, which lies at the upper end of the measured range of Efron et al.59 

In this paper, we investigate the dynamics of tear film and the temperature in a thick rectangular substrate posterior to it. We extend Li and Braun's model16 to include one moving end corresponding to the superior eyelid. The moving end will repeat the blink motion to approximate multiple blink cycles. We investigate the dynamics of the temperature field during and just after a blink, for which we are aware of no experimental data. These dynamics may have an impact on lipid layer dynamics, for example, which have temperature dependent properties.5,45 From this model, we will compute the thermal dynamics posterior to the tear film in the transients around blinking. Previous heat transfer models were for steady states and experimental results do not give temperatures within the cornea and aqueous humor; thus, our results will shed new light on this process. We shall also study tear film dynamics (for the aqueous layer) during the blink cycle. Previous work studied when periodic solutions could occur in the absence of heat transfer;31,32 this periodicity was hypothesized to indicate when a partial blink was fluid-dynamically equivalent to a complete blink. This previous work found that for blinks that were sufficiently close to, but not completely closed, the fluid dynamics did become periodic; this norm of the difference between successive film profiles decreased dramatically at a critical value of blink amplitudes leaving a range of amplitudes near a full blink where periodic solutions occurred. In our current model, we find a slightly smaller range of blink amplitudes where the solutions are periodic. It is beyond the scope of this paper to include lipid layer dynamics,36 but our current work is a significant step in that direction. For both the temperature and tear film dynamics, we make direct comparison with in vivo experimental results where possible. We will also discuss some criteria for partial blinks to maintain the periodic film dynamics subject to the evaporative loss.

We give a brief derivation of the model and introduce the numerical methods in Sec. II, results will be discussed in Sec. III, and compared with experimental data, followed by discussion and conclusion in Secs. IV and V.

A sketch of the mathematical model for the PCTF and a rectangle posterior to it representing the underlying eye is shown in Figure 1. The coordinate direction (x, y) and velocity components (u′, v′) are along and perpendicular to the corneal surface. We assume that the cornea is flat, following many studies (e.g., Berger and Corrsin;23 detailed justification is given in Braun et al.18) Primes indicate dimensional variables.

FIG. 1.

A sketch of pre-corneal tear film with heat transfer from the cornea indicating important mathematical quantities. The dimensional upper lid location is X(t); this end moves in blink cycle models. The grey rectangle represents a pre-existing fluid layer under the lid with thickness

$h^{\prime }_e$
he⁠.

FIG. 1.

A sketch of pre-corneal tear film with heat transfer from the cornea indicating important mathematical quantities. The dimensional upper lid location is X(t); this end moves in blink cycle models. The grey rectangle represents a pre-existing fluid layer under the lid with thickness

$h^{\prime }_e$
he⁠.

Close modal

The tear film/cornea interface is located along the y = 0 edge of the rectangle representing the interface of the tear film and the cornea. The tear film domain varies during a blink and is given by X(t) ⩽ xL as shown in Figure 1, with film thickness given by y = h(x, t). The origin is located in the center of the open eye domain. The function X(t) specifies the upper lid position. When the interblink period is studied, X = −L for the duration of the interblink. The rectangular region bounded by −L < x < L and

$y^{\prime }_c < y^{\prime } < 0$
yc<y<0 is the substrate domain for heat diffusion under the tear film. In this paper, we take
$y^{\prime }_c \gg d^{\prime }$
ycd
, where d represents the typical thickness of the tear film, making the substrate thick compared with the tear film.16,54 In physical terms,
$y^{\prime }_c$
yc
is chosen large enough so that it corresponds to a location inside the aqueous humor. The tear film is modeled as a single layer of incompressible Newtonian fluid with density ρ, dynamic viscosity μ, specific heat cp and thermal conductivity k. Dimensional parameter values are given in Table I.

Table I.

Dimensional parameter definitions and values. Unless otherwise noted, these values were used to generate the computational results. K and

$\bar{K}$
K¯ are chosen to satisfy the experimentally measured thinning rate (4 μm/min73) of the tear film thickness due to King-Smith et al. α and A* are recovered from the nondimensional parameters in Winter et al.27 

Dimensional parameters
ParameterDescriptionValue
d Characteristic thickness of tear film 5 × 10−6
L Half length of the palpebral fissure 5 × 10−3
$h^{\prime }_e$
he
 
Thickness of tear fluid under eyelid 0.25d to d 
$L^{\prime }_c$
Lc
 
Corneal thickness 5 × 10−4
$y^{\prime }_c$
yc
 
Model depth posterior to film 
$5L^{\prime }_c$
5Lc
 
μ Viscosity42  1.3 × 10−3 Pa s 
σ Surface tension47  0.045 N m−1 
k Tear film thermal conductivity63 (water) 0.68 W m−1 K−1 
cp Specific heat (water) 4184 J kg−1 (C°)−1 
kc Corneal thermal conductivity63  0.58 W m−1 K−1 
ρ Density (water) 103 kg m−3 
${\cal L}_m$
Lm
 
Latent heat of vaporization (water) 2.3 × 106 J kg−1 
$T^{\prime }_s$
Ts
 
Saturation temperature (assumed) 27 °C 
$T^{\prime }_{\text{B}}$
TB
 
Body temperature 37 °C 
g Gravitational acceleration 9.81 m s−1 
κc Thermal diffusivity of cornea63  1.9 × 10−7 m2 s−1 
AHamaker constant27  3.5 × 10−19 Pa m3 
U Representative interblink speed 
$U^{\prime } = 5\ \textrm {mm\,s}^{-1}$
U=5 mm s1
 
α Pressure coefficient for evaporation27  3.6 × 10−2 K Pa−1 
K Non-equilibrium coefficient16  1.5 × 105 K m2 s kg−1 
Dimensional parameters
ParameterDescriptionValue
d Characteristic thickness of tear film 5 × 10−6
L Half length of the palpebral fissure 5 × 10−3
$h^{\prime }_e$
he
 
Thickness of tear fluid under eyelid 0.25d to d 
$L^{\prime }_c$
Lc
 
Corneal thickness 5 × 10−4
$y^{\prime }_c$
yc
 
Model depth posterior to film 
$5L^{\prime }_c$
5Lc
 
μ Viscosity42  1.3 × 10−3 Pa s 
σ Surface tension47  0.045 N m−1 
k Tear film thermal conductivity63 (water) 0.68 W m−1 K−1 
cp Specific heat (water) 4184 J kg−1 (C°)−1 
kc Corneal thermal conductivity63  0.58 W m−1 K−1 
ρ Density (water) 103 kg m−3 
${\cal L}_m$
Lm
 
Latent heat of vaporization (water) 2.3 × 106 J kg−1 
$T^{\prime }_s$
Ts
 
Saturation temperature (assumed) 27 °C 
$T^{\prime }_{\text{B}}$
TB
 
Body temperature 37 °C 
g Gravitational acceleration 9.81 m s−1 
κc Thermal diffusivity of cornea63  1.9 × 10−7 m2 s−1 
AHamaker constant27  3.5 × 10−19 Pa m3 
U Representative interblink speed 
$U^{\prime } = 5\ \textrm {mm\,s}^{-1}$
U=5 mm s1
 
α Pressure coefficient for evaporation27  3.6 × 10−2 K Pa−1 
K Non-equilibrium coefficient16  1.5 × 105 K m2 s kg−1 

The scalings we will use to derive our models are as follows: L is the half width of the palpebral fissure and is used as the length scale in x direction; the characteristic thickness of the tear film away from the ends is d and is used in the positive y direction. Since the cornea thickness is

$L^{\prime }_c$
Lc⁠, it is used as the length scale in the negative y direction. The ratio of the length scales for the film, ε = d/L = 10−3, is the small parameter for lubrication theory. The velocity scale along the film is a representative interblink speed, U = 5 mm/s. This is a representative speed of the film surface that occurs at the end of a blink.6,48 At subsequent times, the speeds are usually lower, but this is a convenient choice for our work. We note that εU is the characteristic speed across the film. The time is scaled by L/U, the pressure p is made non-dimensional with viscous scale μU/(Lε2), the evaporation rate J is non-dimensionalized by scale
$\frac{k}{d^{\prime } {\cal L}_m}(T^{\prime }_{\text{B}}-T^{\prime }_s)$
kdLm(TBTs)
and the temperature T is made non-dimensionalized as
$T = \frac{T^{\prime }-T^{\prime }_s}{T^{\prime }_{\text{B}}-T^{\prime }_s}$
T=TTsTBTs
. Here
$T^{\prime }_B$
TB
is body temperature,
$T^{\prime }_s$
Ts
is the saturation temperature corresponds to saturation pressure at which a tear fluid equilibrates with its vapor phase, k is the thermal conductivity, and
${\cal L}_m$
Lm
represents the latent heat of vaporization per unit mass. The numerical values of the parameters are given in Table II.

Table II.

Non-dimensional parameter definitions and values.

 Dimensionless parameters 
ε 
$\frac{d^{\prime }}{L^{\prime }}$
dL
 
10−3 
Re 
$\frac{\rho U^{\prime } L^{\prime }}{\mu }$
ρULμ
 
20 
Pr 
$\frac{\mu c_P}{k}$
μcPk
 
8.0 
E 
$\frac{k(T^{\prime }_{\text{B}}-T^{\prime }_s)}{{d^{\prime } {\cal L}_m} \epsilon \rho U^{\prime }}$
k(TBTs)dLmερU
 
118.3 
S 
$\frac{\sigma \epsilon ^3}{\mu U^{\prime }}$
σε3μU
 
6.92 × 10−6 
$\bar{K}$
K¯
 
$\frac{kK}{{d^{\prime } {\cal L}_m}}$
kKdLm
 
8.9 × 103 
he 
$\frac{h^{\prime }_e}{d^{\prime }}$
hed
 
δ 
$\frac{\alpha \mu U^{\prime }}{\epsilon ^2L(T^{\prime }_{\text{B}}-T^{\prime }_s)}$
αμUε2L(TBTs)
 
4.66 
A 
$\frac{A^*}{{L^{\prime } d^{\prime }} \mu U^{\prime }}$
A*LdμU
 
2.14 × 10−6 
$\tilde{k}$
k̃
 
k/kc 1.17 
PT 
$\frac{L^{\prime } \kappa _c}{U^{\prime } {{L^{\prime }_c}^2}}$
LκcULc2
 
0.76 
Bi 
$\frac{h_cd^{\prime }}{k}$
hcdk
 
0-0.1 
 Dimensionless parameters 
ε 
$\frac{d^{\prime }}{L^{\prime }}$
dL
 
10−3 
Re 
$\frac{\rho U^{\prime } L^{\prime }}{\mu }$
ρULμ
 
20 
Pr 
$\frac{\mu c_P}{k}$
μcPk
 
8.0 
E 
$\frac{k(T^{\prime }_{\text{B}}-T^{\prime }_s)}{{d^{\prime } {\cal L}_m} \epsilon \rho U^{\prime }}$
k(TBTs)dLmερU
 
118.3 
S 
$\frac{\sigma \epsilon ^3}{\mu U^{\prime }}$
σε3μU
 
6.92 × 10−6 
$\bar{K}$
K¯
 
$\frac{kK}{{d^{\prime } {\cal L}_m}}$
kKdLm
 
8.9 × 103 
he 
$\frac{h^{\prime }_e}{d^{\prime }}$
hed
 
δ 
$\frac{\alpha \mu U^{\prime }}{\epsilon ^2L(T^{\prime }_{\text{B}}-T^{\prime }_s)}$
αμUε2L(TBTs)
 
4.66 
A 
$\frac{A^*}{{L^{\prime } d^{\prime }} \mu U^{\prime }}$
A*LdμU
 
2.14 × 10−6 
$\tilde{k}$
k̃
 
k/kc 1.17 
PT 
$\frac{L^{\prime } \kappa _c}{U^{\prime } {{L^{\prime }_c}^2}}$
LκcULc2
 
0.76 
Bi 
$\frac{h_cd^{\prime }}{k}$
hcdk
 
0-0.1 

1. General problem description

After scaling as above, we obtain dimensionless (unprimed) variables and equations. For the flow in the film on X(t) < x < 1 and 0 ⩽ yh(x, t), we have conservation of mass, momentum (in the x and y directions) and energy, respectively:

\begin{eqnarray}& u_x + v_y = 0 ,\end{eqnarray}
ux+vy=0,
(1)
\begin{eqnarray}& \epsilon ^2\text{Re}(u_t + uu_x+vu_y) = -p_x + \epsilon ^2u_{xx} + u_{yy} ,\end{eqnarray}
ε2Re(ut+uux+vuy)=px+ε2uxx+uyy,
(2)
\begin{eqnarray}& \epsilon ^4\text{Re}(v_t + uv_x+vv_y) = -p_y + \epsilon ^4v_{xx} + \epsilon ^2v_{yy} ,\end{eqnarray}
ε4Re(vt+uvx+vvy)=py+ε4vxx+ε2vyy,
(3)
\begin{eqnarray}& { \epsilon ^2 \text{RePr}(T_t + u T_x + v T_y ) } = \epsilon ^2 T_{xx} + T_{yy}.\end{eqnarray}
ε2RePr(Tt+uTx+vTy)=ε2Txx+Tyy.
(4)

The Reynolds number Re = ρUL′/μ ≈ 20, and together with ε = 10−3, we have the reduced Reynolds number ε2Re ≪ 1. For the energy equation, the Prandtl number Pr = μc/k ≈ 5 at the ocular surface temperature, where c is the specific heat and we have assumed water properties for the aqueous layer of the tear film. The quantity ε2RePr is also small, so we may safely neglect the inertial and convective terms in the model derivation below.

We assume that the bottom surface posterior to the tear film is both impermeable and no slip so that u = v = 0 there. On the free surface, standard conditions apply: the kinematic condition, and normal and tangential stress conditions;74,75 a constitutive equation for evaporation that is linear the film pressure and temperature at the free surface;27,51 and the transport of an insoluble surfactant with a linearized constitutive equation for the surface tension.31,75 The normal stress contribution includes a contribution from van der Waals forces that model the wettable corneal surface as described previously.17,27

In this work, we assume that the lipid layer acts like an insoluble surfactant.24–26,35 Its effect on the free surface may be further simplified to the USL corresponding to a strong Marangoni effect;24,31,32 a detailed derivation of the USL has been given in Ref. 31 and we give the resulting condition on the surface velocity below. We assume that the evaporation rate is unaffected by the presence of the insoluble surfactant.

Under the moving eyelid, corresponding to the interval −Lx < X(t), we assume that there is a thin fluid layer of constant thickness specified by

$h^{\prime }_e=h_e d^{\prime }$
he=hed⁠.24,32,35 Here he is dimensionless, and indicates the relative thickness of the gap compared with the characteristic film thickness. In this gap, a Couette flow occurs between the moving lid and the stationary corneal surface. Both the eyelid and ocular surface on either side of the gap have been idealized to be flat, smooth, impermeable, and no slip surfaces. Regarding heat transfer, we assume that heat is conducted across the thin fluid gap with corneal surface temperature T(x, 0, t) on the corneal side and (constant) body temperature
$T^{\prime }_B$
TB
on the eyelid surface.

At the film ends, we need boundary conditions as well. At those ends, the tear film is assumed to be pinned at a thickness h0d; h0 is a dimensionless measure of the thickness at which the tear film is pinned to the eyelid relative to the characteristic film thickness. In the interval hed < y < h0d, the tear film velocity components are equal to the motion of the ends. In the interval 0 < y < hed, we may specify the velocity components; however, only the depth-average of a such a velocity profile (or the equivalent flux) will be applied to the resulting thin film equations. We hypothesize two components of the flux. One component arises from the Couette flow in the gap between the moving (superior) eyelid and the corneal surface, which is proportional to the eyelid (domain end) speed.24–26,32,35 The other is from the lacrimal gland supply and punctal drainage, which is not proportional to the end motion, but its timing should be chosen properly in relation to the eyelid motion.32,33 We specify the end fluxes and their components below.

2. Thin film approximation

After dropping terms containing ε in Eqs. (1)–(3), we have the following leading-order problem:

\begin{equation}u_x+v_y = 0, \ \ u_{yy} -p_x = 0, \ \ p_y = 0.\end{equation}
ux+vy=0,uyypx=0,py=0.
(5)

At the free surface of the film y = h(x, t), we have the kinematic and normal stress conditions, respectively,

\begin{equation}h_t + u h_x + E J = v,\end{equation}
ht+uhx+EJ=v,
(6)
\begin{equation}p = -S h_{xx} - A/h^3 .\end{equation}
p=ShxxA/h3.
(7)

Here E gives the size of the nondimensional evaporative mass flux J from the film surface and A is the nondimensional Hamaker constant in the unretarded van der Waals force. The reduced surface tension number S measures the effect of surface tension or capillary force. Thus the model includes the evaporative effect from the aqueous layer and the van der Waals forces52,53 representing the wettable corneal surface. In the USL, the tangential stress condition at the free surface is replaced with a boundary condition on the surface velocity that may be interpreted as the linear stretching of the surface:

\begin{equation}u_s \equiv u(x,h,t) = \frac{1-x}{1-X} X_t ,\end{equation}
usu(x,h,t)=1x1XXt,
(8)

The thermal energy balance and the constitutive equation for the evaporative mass flux are given respectively by

\begin{equation}J +T_y + \text{Bi}(T-T_{\infty }) = 0 ,\end{equation}
J+Ty+Bi(TT)=0,
(9)
\begin{equation}\bar{K}J = \delta p + T .\end{equation}
K¯J=δp+T.
(10)

Here Bi is the Biot number indicating the ratio of the heat transfer resistance from convection outside the film to conduction across the tear film.

$\bar{K}$
K¯ and δ are nondimensional parameters relating the evaporation rate and the temperature and pressure at the free surface.

At leading order, the energy equation (4) reduces to

\begin{equation}T_{yy} = 0.\end{equation}
Tyy=0.
(11)

Solving for the linear temperature field in terms of the nondimensional evaporative mass flux J yields

\begin{equation}T(x,y,t) = -\frac{J+\text{Bi}(T(x,0,t)-T_{\infty })}{1+\text{Bi} h}y+T(x,0,t), \ \ -1 < x < 1 .\end{equation}
T(x,y,t)=J+Bi(T(x,0,t)T)1+Bihy+T(x,0,t),1<x<1.
(12)

J is found by evaluating T at y = h using Eq. (12), substituting into Eq. (10) and solving for J. The governing equation for the film thickness may then be derived in the usual way from the leading order problem, yielding16 

\begin{equation}h_t + E J + q_x = 0 , \ {X(t) < x < 1} ,\ t> 0 ,\end{equation}
ht+EJ+qx=0,X(t)<x<1,t>0,
(13)

where

\begin{equation}q = \frac{h^3}{12} (S h_{xx} + Ah^{-3})_x + X_t \frac{1-x}{1-X} \frac{h}{2},\end{equation}
q=h312(Shxx+Ah3)x+Xt1x1Xh2,
(14)

and

\begin{equation}J = \left[\frac{T_0 + \text{Bi} T_{\infty } h}{1+ \text{Bi} h} - \delta (S h_{xx} +Ah^{-3})\right]\left(\bar{K} + \frac{h}{1+\text{Bi} h}\right)^{-1}.\end{equation}
J=T0+BiTh1+Bihδ(Shxx+Ah3)K¯+h1+Bih1.
(15)

Here

\begin{equation}T_0=T(x,0,t)=\tilde{T}(x,0.t)\end{equation}
T0=T(x,0,t)=T̃(x,0.t)
(16)

is the temperature at the tear film/cornea interface and is a dependent variable to be found as part of the solution.

In the rectangular region posterior to the tear film, the temperature changes are only due to heat diffusion. The small effect of fluid motion inside the eye is negligible.67 Let

$\tilde{T}(x,y,t)$
T̃(x,y,t) denote the nondimensional temperature posterior to the film
$\tilde{y_c} \le \tilde{y} \le 0, \ -1 \le x\le 1$
yc̃ỹ0,1x1
, subject to

\begin{equation}\tilde{T}_t = P_T \left[ \left(\frac{L^{\prime }_c}{L^{\prime }}\right)^2 \tilde{T}_{xx} + \tilde{T}_{yy} \right].\end{equation}
T̃t=PTLcL2T̃xx+T̃yy.
(17)

The heat domain and film domain are connected through a time-dependent boundary y = 0, X(t) ⩽ xL by applying continuity of temperature and heat flux. The boundary conditions in their dimensionless forms are

\begin{equation}u = v = 0, \ \ T_0 = \tilde{T}, \ \ \text{where} \ \ T_0 = T(x,0,t),\end{equation}
u=v=0,T0=T̃,whereT0=T(x,0,t),
(18)
\begin{equation}\tilde{T}_{\tilde{y}} + \frac{\tilde{k}L^{\prime }_c}{d^{\prime }} \frac{J + \text{Bi}( T_0 - T_{\infty })}{1+\text{Bi} h} = 0, \ \ X(t) < x < 1,\end{equation}
T̃ỹ+k̃LcdJ+Bi(T0T)1+Bih=0,X(t)<x<1,
(19)
\begin{equation}\tilde{T}_{\tilde{y}} - \frac{\tilde{k}L^{\prime }_c}{h_e d^{\prime }} (T_{\text{B}} - T_0) = 0, \ \ -1 < x < X(t).\end{equation}
T̃ỹk̃Lched(TBT0)=0,1<x<X(t).
(20)

Depending on whether a particular location in x is covered by the upper eyelid or not, we have two different heat transfer rates. As sketched in Figure 2, the two boundary conditions can be joined together smoothly with

\begin{equation}\psi (x) = \frac{1}{2} \left[ 1 + \tanh {\left(\frac{x-X(t)}{x_w}\right)}\right],\end{equation}
ψ(x)=121+tanhxX(t)xw,
(21)

Our choice of ψ(x) is motivated by the need to smoothly join the two different boundary conditions (19)-(20) on y = h, so that we can resolve the computations with increasing grid refinement. The hyperbolic tangent function serves this purpose well. Thus, on −1 ⩽ x ⩽ 1, we use

\begin{equation}\tilde{T}_{\tilde{y}}(x,0,t) + (1-\psi ) \frac{L^{\prime }_c \tilde{k}}{h^{\prime }_e d^{\prime }} (T_0 - T_{\text{B}} ) + \psi \text{Bi} (T_0-T_{\infty }) = 0.\end{equation}
T̃ỹ(x,0,t)+(1ψ)Lck̃hed(T0TB)+ψBi(T0T)=0.
(22)

Here xw = 0.01 is chosen to be small, so that there is a narrow band for the transition between the two conditions.

FIG. 2.

Neumann conditions corresponding to open air and upper eyelid can be joined smoothly with a hyperbolic tangent function. The sketch shown in not to scale; the width of the transition from one boundary condition is only about 1% of the domain length.

FIG. 2.

Neumann conditions corresponding to open air and upper eyelid can be joined smoothly with a hyperbolic tangent function. The sketch shown in not to scale; the width of the transition from one boundary condition is only about 1% of the domain length.

Close modal

The other three edges of the rectangle are prescribed to be at body temperature, with boundary conditions

$T^{\prime }_{B} = 37\,^{\circ }\text{C}$
TB=37°C⁠; after scaling, we have

\begin{equation}\tilde{T}(x,\tilde{y_c},t) = \tilde{T}(\pm 1,\tilde{y},t ) = 1.\end{equation}
T̃(x,yc̃,t)=T̃(±1,ỹ,t)=1.
(23)

This assumption is based on the fact that blood flow in the lids and the sclera all around the eye help to heat the deeper interior of the eye, and thus maintain a relatively constant temperature that is close to body temperature.

The eyelid motion of a blink is solved on a time-periodic domain with specified film thickness and flux at each end. It is assumed that only the left end X(t) of the domain moves; the motion corresponds to the superior eyelid movement over blink cycles. A mathematical fit to the lid position at the center of palpebral fissure was first obtained by Berke and Mueller.76 Other lid motion functions have been developed,24,26,32 and we use those from Heryudono et al.32 They modified the function used by Berke and Meuller76 to describe both partial and full blinks; that lid motion function, is referred to as realistic lid motion here. The lid motion is defined on the interval −1 ⩽ X(t) ⩽ 1 − 2λ and the tear film is on the interval X(t) ⩽ x ⩽ 1. Here λ is the fraction of the fully open region that is still open at the end of the downstroke, when the domain is at its smallest. Using our non-dimensionalization, the duration of the upstroke is Δtco = 0.176; the downstroke, Δtoc = 0.082, and the interblink time, Δto = 5. Thus, the blink cycle duration to be used in our experiment representing a normal eye blink is Δtbc = 5.258. A detailed expression of X(t) is given in Appendix  A.

Braun and King-Smith31 explored full blink cycles using an ideal case: the sinusoidal lid motion

\begin{equation}X(t) = (1-\lambda )\cos {t} - \lambda .\end{equation}
X(t)=(1λ)costλ.
(24)

We will follow the approaches in Refs. 31 and 32 and consider film with λ = 0.1 as being fully closed for this paper, which we refer to as a “full” blink. On the other hand, any blink with λ > 0.1 will be considered as a “partial” blink, as observed from many subjects participating the experiments of Efron et al.59 

The flux boundary conditions of fluid into and out of the tear film domain model the flow along the lids during blink cycles as described by Doane.39 The flux conditions are

\begin{equation}q(X(t),t) = X_t h_0 + Q_{\text{top}}, \ \ q(1,t) = -Q_{\text{bot}} ,\end{equation}
q(X(t),t)=Xth0+Qtop,q(1,t)=Qbot,
(25)

where Qtop and Qbot represent specified fluxes into the domain at the moving and stationary ends respectively. Jones et al.24 observed that a significant volume of tear fluid is stored under the superior and inferior eyelids, and thus as the superior lid moves up or down, there is a flux of fluid from under the moving superior eyelid. They assumed that the layer thickness is a constant

$h^{\prime }_e$
he⁠, and hence Qtop relative to the lid is proportional to the eyelid velocity. We refer to this flux condition as the “Flux Proportional to Lid Motion” (FPLM) condition, and non-dimensionally it takes the form

\begin{equation}Q_{\text{top}} = -X_t h_e/2 , \ \ Q_{\text{bot}} = 0.\end{equation}
Qtop=Xthe/2,Qbot=0.
(26)

Heryudono et al.32 explored the film dynamics’ dependence on he. We note that if Qtop = 0, that the volume will be constant so that (25) become no flux boundary conditions. By analogy, the correct condition for exposing film a film of constant thickness from under the moving lid may be shown to be (26) because of the exposed tear volume.24,32 Unless otherwise specified, we use he = 1 in our computations, corresponding to

$h^{\prime }_e=5\,\,\mu$
he=5μm. An example of the FPLM boundary conditions is shown in Figure 3.

FIG. 3.

An illustration of realistic lid motion X(t), lid speed Xt(t) and boundary fluxes Qtop and Qbot for a full blink with FPLM boundary condition. The solid curve represents the lid position X(t) during the blink, which starts and ends slightly off x = 1 (lower eyelid), meaning a fraction λ of the domain remains open to air over multiple blinks. The dashed line gives the corresponding lid velocity Xt(t). The dotted line gives the influx or efflux of fluid from the superior eyelid Qtop due to our assumption on a preexisting fluid layer between the lid and cornea as specified by Eq. (26). The dash-dot line represents the flux of fluid from the lower eyelid Qbot, which will remain zero for the FPLM condition, but will be positive as a result of lacrimal gland supply in the FPLM+ condition later. Note that the t axis is not to scale.

FIG. 3.

An illustration of realistic lid motion X(t), lid speed Xt(t) and boundary fluxes Qtop and Qbot for a full blink with FPLM boundary condition. The solid curve represents the lid position X(t) during the blink, which starts and ends slightly off x = 1 (lower eyelid), meaning a fraction λ of the domain remains open to air over multiple blinks. The dashed line gives the corresponding lid velocity Xt(t). The dotted line gives the influx or efflux of fluid from the superior eyelid Qtop due to our assumption on a preexisting fluid layer between the lid and cornea as specified by Eq. (26). The dash-dot line represents the flux of fluid from the lower eyelid Qbot, which will remain zero for the FPLM condition, but will be positive as a result of lacrimal gland supply in the FPLM+ condition later. Note that the t axis is not to scale.

Close modal

To describe the lacrimal gland supply of new tear fluid and the drainage of excess tears, Heryudono et al.32 modified the FPLM condition by adding new terms representing influx from the lacrimal gland and outflux through the puncta. The specific functions are given in Appendix  A. We refer to this flux condition as the FPLM+ condition.32 

The problem to be solved is for the tear film thickness given by Eqs. (13)–(15) the temperature inside the model eye domain given by Eqs. (17), (22), and (23) and the continuity of temperature at y = 0, Eq. (16). We impose the following boundary conditions:

\begin{equation}h(1,t) = h_0, \ \ h(X(t),t) = h_0,\end{equation}
h(1,t)=h0,h(X(t),t)=h0,
(27)

and

\begin{equation}q(1,t) = - Q_{\text{bot}}, \ \ q(X(t),t) = X_t h_0 + Q_{\text{top}} ,\end{equation}
q(1,t)=Qbot,q(X(t),t)=Xth0+Qtop,
(28)

where Qbot and Qtop depend on the choice of flux conditions. The initial condition for the film is specified as

\begin{equation}h(x,0) = h_{\text{min}} + (h_0 - h_{\text{min}}) \left[1- 2\frac{1-x}{1-X(0)}\right]^m .\end{equation}
h(x,0)=hmin+(h0hmin)121x1X(0)m.
(29)

Here m characterizes the shape of the polynomial and hmin is computed such that the volume V0 = V(0) under the curve is a specified value. The temperature field is initialized with body temperature:

\begin{equation}\tilde{T}(x,\tilde{y},0) = 1.\end{equation}
T̃(x,ỹ,0)=1.
(30)

We used the method of lines with a Chebyshev spectral collocation discretization in space to achieve high accuracy.77 We found better performance in avoiding the fourth-order derivative in space that appears in the thin film equation by introducing the pressure p(x, t) = −(Shxx + Ah−3) as an additional dependent variable. This reduces the highest order derivative to two, but doubles the number of unknowns to solve for in the thin film. The spatial discretization creates a system of differential-algebraic equations (DAE) in time. Details of the method are given in Appendix  B.

We begin with results using sinusoidal lid motion (24), where the FPLM condition will be used. We then proceed to realistic lid motion (A1), where both the FPLM condition and FPLM+ condition will be considered.

In this section, we will explore the dynamics resulting from sinusoidal lid motion, given by (24) and used by Braun and King-Smith,31 in a more general model including evaporation and heat transfer posterior to the PCTF. The parameters used in this section are S = 7 × 10−6, V0 = 1.576, λ = 0.1, m = 4, and h0 = 13. Note that λ = 0.1 represents 10% of domain left open at the end of a blink cycle; thus, we are dealing with full blinks in this section, unless otherwise specified. To best simulate the experimental results in our computations, we use Bi = 0.0009 after Li and Braun;16 in later sections we vary it and will discuss its effect.

We employ the FPLM flux condition (26) to model a flux entering or leaving the film due to the covering or uncovering of the preexisting fluid layer under the moving lid.24 Figure 4 shows the film profile for a blink cycle with he = 1. During the upstroke, the thickness of the tear film is deposited at more than a plausible minimum thickness because of a fluid influx from the FPLM boundary condition. This is in contrast to the no flux boundary condition, where the deposited tear film may be thin enough to cause film rupture (also known as break up; not shown).24 In Figure 5, we show the film thickness profiles during the interblink (open phase) over four successive blink cycles. After the first blink cycle, the global minimum is near the moving superior lid in accord with many other works.24,26,32,34,35 After subsequent blink cycles, the global minimum is located near the stationary inferior lid. We believe that this is due to the sinusoidal lid motion, and adding evaporation to this simplified motion reveals previously unseen limitations from it. Turning to the total volume of tear fluid, previous work has found that in the absence of evaporation, the total volume (area under h) would be equal each after each period of t = 5.28 time units to within numerical error (not shown).31,32 When the tear film evaporates, some volume is lost during each cycle; it can be seen that most of the evaporative loss occurs from the inferior (stationary) end for sinusoidal lid motion. The total volume of the tear fluid is plotted in Figure 6. The amount of fluid oscillates because of the FPLM boundary condition adding or removing fluid depending on the sign of Xt, but on the average the volume decreases over time due to evaporation.

FIG. 4.

Film thickness h(x, t) at t = nπ/4 for n = 0, 1, ⋅ ⋅ ⋅, 8 with λ = 0.1, V0 = 1.576, and FPLM boundary conditions.

FIG. 4.

Film thickness h(x, t) at t = nπ/4 for n = 0, 1, ⋅ ⋅ ⋅, 8 with λ = 0.1, V0 = 1.576, and FPLM boundary conditions.

Close modal
FIG. 6.

Tear film total volume decreases for four successive periods; circles are at the same time offset within one period for all four periods.

FIG. 6.

Tear film total volume decreases for four successive periods; circles are at the same time offset within one period for all four periods.

Close modal
FIG. 5.

Film profiles shown for four successive periods when the eye is completely open, subject to the FPLM condition.

FIG. 5.

Film profiles shown for four successive periods when the eye is completely open, subject to the FPLM condition.

Close modal

The temperature field corresponding to the first two blink cycles is shown in Figure 7. The correlation between the moving end of the film and temperature is clear. During the upstroke, the open eye is cooled by the air outside. During the downstroke, the superior eyelid supplies heat that rapidly warms up the eye close to the tear film, while leaving a visible cool “island” inside the eye. We also see that a very similar temperature profile is established inside the eye at the end of each blink. Subsequent blink cycles may approach periodic temperature behavior if the blinks are sufficiently full.

FIG. 7.

Temperature profile for a complete blink cycle with λ = 0.1, V0 = 1.576, sinusoidal lid motion, and FPLM boundary condition. The vertical sides and bottom edge are always warmest. To our knowledge there is no available experimental data measuring the temperature inside the cornea. Here our simulation is meant to suggest qualitative behavior. We defer a quantitative comparison with measured GCC temperature later for the real blink case.

FIG. 7.

Temperature profile for a complete blink cycle with λ = 0.1, V0 = 1.576, sinusoidal lid motion, and FPLM boundary condition. The vertical sides and bottom edge are always warmest. To our knowledge there is no available experimental data measuring the temperature inside the cornea. Here our simulation is meant to suggest qualitative behavior. We defer a quantitative comparison with measured GCC temperature later for the real blink case.

Close modal

According to Doane,39 human eyes usually have several partial blinks between a pair of full blinks. For partial blinks, the superior eyelid stops before reaching the inferior eyelid at the end of the downstroke. This behavior is modeled by setting λ > 0.1 in our model. Braun and King-Smith31 first found that periodicity in film thickness could be achieved for partial blinks over multiple blink cycles. They explored the sensitivity of periodic tear film dynamics to different λ for sinusoidal lid motion (24) and found that there exists a critical value for λ, somewhere in the range

$\left[\frac{1}{8},\frac{1}{6}\right]$
18,16⁠, below which a periodic solution could be achieved. They attributed this behavior to the fact that tear fluid is a slow viscous flow, and the capillary-driven film shape formed near the stationary end during each blink could be eliminated at the end of downstroke. Heryudono et al.32 explored this phenomenon further with realistic lid motion (A1), and found a similar conclusion. We refer to the phenomenon that certain partial blinks could achieve periodic tear film dynamics as full blinks as “full blink equivalence.” The aforementioned work did not consider evaporation, which can play a significant role in film deformation,28 especially when the blink rate is slow. To explore the full blink equivalence for the model including evaporation, we need to compensate for the evaporation loss by supplying the right amount of fluid through the superior eyelid. We then test the sensitivity of periodic solutions to different values of λ by measuring the difference between film profiles from successive blinks in the Euclidean norm for successive blink cycles.

In Figure 6, the volume decreases at an average rate that is effectively constant over several blink cycles. The circles in the figure indicate volumes separated by one period, and the dashed line illustrates the average rate of volume decrease. We compute a compensation flux Qc needed to replace the evaporative loss based on the average rate of decrease of V for the FPLM case. We then modify the flux through the superior eyelid to Qtop = −Xthe/2 + Qc to be used in Eq. (28); we refer to this approach as the compensated FPLM (cFPLM) condition. With this flux condition, we guarantee that the volume V under the film remains periodic over multiple blink cycles. To validate this argument, Figure 8 shows the effectiveness of the cFPLM condition: the volume difference over one period, given by |V(α + 2kπ) − V(α)|, indicates that up to four digits of precision can be obtained with this method of evaporation compensation for several blink cycles.

FIG. 8.

Volume difference between subsequent blinks after cFPLM evaporation compensation at the different time levels. Note the small scale for the ordinate.

FIG. 8.

Volume difference between subsequent blinks after cFPLM evaporation compensation at the different time levels. Note the small scale for the ordinate.

Close modal

Figures 9(a) and 9(b) show the film thickness profiles during successive interblinks (domain is fully open) over four consecutive blink cycles with λ = 0.125 and λ = 0.15, respectively. Both represent partial blinks according to our definition. We observe that the tear film thickness is periodic for λ = 0.125. For a slightly larger λ = 0.15, the solution has a thin region near the right meniscus, the so-called black line, that deepens after each blink cycle. In this latter case, the eyelid does not push down enough to sweep away the deformation as a combined result of surface tension and evaporation during each blink, so the deformation remains and accumulates over each blink.

FIG. 9.

In (a)–(c), film thickness at the start of successive interblinks are shown for different values of λ. Thickness profiles are shown after four successive periods, specifically, t = π, 3π, 5π, 7π using the FPLM condition. We see periodic solutions at λ = 0.125 and non-periodic solutions at λ = 0.15. Different film behaviors are shown for λ = 0.15 with no evaporation or evaporation compensation. In (d), full blink equivalence is seen when the norm of the difference between successive film thickness profiles (Δ, see text) decreases suddenly (right to left).

FIG. 9.

In (a)–(c), film thickness at the start of successive interblinks are shown for different values of λ. Thickness profiles are shown after four successive periods, specifically, t = π, 3π, 5π, 7π using the FPLM condition. We see periodic solutions at λ = 0.125 and non-periodic solutions at λ = 0.15. Different film behaviors are shown for λ = 0.15 with no evaporation or evaporation compensation. In (d), full blink equivalence is seen when the norm of the difference between successive film thickness profiles (Δ, see text) decreases suddenly (right to left).

Close modal

The average film volume V is kept constant over successive blink cycles in order to facilitate comparison the case without evaporation.31 Since the evaporation takes place all over the film, and the compensation flux Qc enters only through the moving end (superior eyelid), we expect different film dynamics compared with those without evaporation. We computed the film profiles over four consecutive blink cycles for these two cases with λ = 0.15, as shown in Figures 9(b) and 9(c). The relatively small opening of λ = 0.15, just 15% of fully open, does not give a periodic solution if one simply replaces the evaporative loss; the domain must close further to eliminate the film thickness variation seen in these two figures. The results as a function of λ are illustrated in Figure 9(d), where we use Δ = ‖h(x, 2tbc + tco) − h(x, tbc + tco)‖2 to measure the deviation between solutions over the second blink cycle. In the presence of evaporation compensation, we find that Δ shifts slightly rightward, leaving a wider λ range in which full blink equivalence is observed. With or without compensation for evaporative water loss, the presence of evaporation of water from the tear film requires that blinks be more complete to achieve periodic solutions compared with previous models without evaporation.32 

In this section, we study the tear film with realistic lid motion given by Eq. (A1) with either the FPLM condition or the FPLM+ condition. In this case, the lid motion of the blink takes up about 5% of the blink cycle, while the remaining time corresponds to the lids being stationary and fully open during the interblink. Heryudono et al.32 observed that periodicity could be obtained if the ends closed sufficiently during the blink cycle. In the following, we will test for full blink equivalence with a model including evaporation. Parameters that will be used in this section are S = 7 × 10−6, Bi = 0.0009, V0 = 2.576, λ = 0.1, m = 4, and h0 = 13, unless otherwise specified. We compute the temperature field corresponding to both the model eye domain and the OST for this realistic lid motion, and will compare the OST with the experimental results of Efron et al.59 We will discuss the film thinning rate and its space dependence at the end of this section.

1. FPLM condition

We start by considering the time sequence of a single blink cycle, where the lid motion is given by Eq. (A1). We divide a complete blink into upstroke, interblink (open) and downstroke phases. Figure 10 depicts film profiles in each phase. During the upstroke (Figure 10(a)) the FPLM condition provides influx/outlfux through the moving end, which serves to avoid film rupture and help the film thickness to distribute uniformly.24,32 During the interblink (Figure 10(b)), both evaporation and surface tension work on the film and influence its shape. The evaporation takes place all over the film, but is most obvious in the middle of the film where there is a nearly uniform drop in film thickness. The downstroke removes some fluid from the tear film and helps refresh the film thickness for the next upstroke. The surface tension tends to increase the minima of the film at both ends, where it competes with evaporation to give two little bumps and then two deep valleys corresponding to the black lines.17,22 From Figure 10(c), it can also be observed that the film deformation from surface tension can be eliminated if the lid closes sufficiently.

FIG. 10.

A complete blink cycle with realistic lid motion as given in Eq. (A1). Most evaporation happens during the interblink.

FIG. 10.

A complete blink cycle with realistic lid motion as given in Eq. (A1). Most evaporation happens during the interblink.

Close modal

We show the correlation between the film and the temperature for a full blink in Figure 11. The temperature field behaves differently from that of sinusoidal lid motion: while sinusoidal lid motion keeps one end always in motion, the realistic lid motion models a normal blink, where the eyelids are fixed fully open for a period of time much longer than that of upstroke and downstroke motion. As a result, most cooling of the eye happens during the open phase. As the lid moves down, heating from the superior eyelid rapidly warms a narrow boundary layer at the anterior side of the model eye domain just underneath the film. The short downstroke phase traps a cool “island” inside the model eye domain as in the sinusoidal case; the island is visible as the isolated interior dark patch (blue online) at t = 5.24. The process may again reach a periodic state over multiple blink cycles.

FIG. 11.

Film and temperature profiles of a blink cycle with realistic lid motion, time levels are selected to reflect the three phases in the blink. In the upstroke, t = 0.12; two time levels in the interblink, t = 2.176 and t = 5.176; in the downstroke, t = 5.24.

FIG. 11.

Film and temperature profiles of a blink cycle with realistic lid motion, time levels are selected to reflect the three phases in the blink. In the upstroke, t = 0.12; two time levels in the interblink, t = 2.176 and t = 5.176; in the downstroke, t = 5.24.

Close modal

Figure 12(a) suggests that the volume of the tear film decreases at an effectively constant rate if we sample the volume at the same point in the blink cycle (this corresponds to the circles in Figure 12(a)). By assuming a constant influx Qc from beneath the superior eyelid, the volume loss from the evaporation is replaced over each blink cycle. With such compensation, we can test for the periodicity of partial blinks using realistic lid motion. Figures 12(b) and 12(c) show that periodic dynamics can be achieved for λ = 0.1 but not for λ = 0.125. Note that λ = 0.1 corresponds to a “full” blink (by definition), so for realistic lid motion we claim that a blink must effectively be complete for periodic tear film dynamics. Figure 12(d) shows the periodicity test as a function of λ for both the sinusoidal lid motion and realistic lid motion. Here we choose Δ = ‖h(x, 3tbc + tco) − h(x, 2tbc + tco)‖2 to measure the deviation between thickness profiles over successive blink cycles. It is clearly seen that the realistic lid motion shows greater sensitivity to λ than the sinusoidal lid motion. We use Δ = 10−2 as a threshold to set the boundary between periodic and non-periodic behaviors. The sinusoidal lid motion, which may be thought of as representing continuous blinking, suggests that periodic behavior can be achieved with λ ⩽ 0.125, so full blink equivalence is possible for some partial blinks. The realistic lid motion, which represents normal spontaneous eye blinking, suggests that only a full blink could achieve periodic dynamics.

FIG. 12.

Top left: Tear film volume with evaporation and FPLM condition. Top right and bottom left: Periodic and non-periodic solutions for different λ with S = 7 × 10−6, V0 = 1.576. Five film profiles are shown at the specific time of t = 2.716 in each of blink cycle. Bottom right: Detection of periodic film behaviour for different values of λ are shown for sinusoidal motion and realistic motion, with S = 7 × 10−6, Bi = 0.0009, V0 = 1.576, m = 4, and h0 = 25.

FIG. 12.

Top left: Tear film volume with evaporation and FPLM condition. Top right and bottom left: Periodic and non-periodic solutions for different λ with S = 7 × 10−6, V0 = 1.576. Five film profiles are shown at the specific time of t = 2.716 in each of blink cycle. Bottom right: Detection of periodic film behaviour for different values of λ are shown for sinusoidal motion and realistic motion, with S = 7 × 10−6, Bi = 0.0009, V0 = 1.576, m = 4, and h0 = 25.

Close modal

2. FPLM+ boundary condition

In this section we will consider the FPLM+ condition (A2)-(A3), where the FPLM condition and Gaussian fluxes from the ends are combined to represent the fluxes that result from lacrimal gland supply and punctal drainage that occur during a blink. Under the FPLM+ condition we will examine full blink equivalence, the surface temperature T(x, h(x, t), t), the GCC temperature T(0, h(0, t), t), and the film thinning rate

$h^{\prime }_{t^{\prime }}(t^{\prime })$
ht(t)⁠. Our simulation results will be compared with the experiment of Efron et al.59 for the GCC temperature as a function of time.

As shown in Figure 13(a), the net effect of the influx due to the lacrimal gland and efflux due to the puncta is to supply fluid into the tear film when the eye is nearly open. However, this flux does not adequately compensate for the evaporation during each cycle; thus we computed the additional flux Qc needed for the compensation as in Sec. III B 1. We call this case the compensated FPLM+, or cFPLM+, boundary condition. Figure 13(b) shows the result of seeking periodicity using cFPLM+ boundary conditions; when Δ = ‖h(x, 2tbc + tco) − h(x, tbc + tco)‖2 drops to values below 10−2, we declare the solution to be periodic. We find that both cFPLM and cFPLM+ conditions yield similar results for the size of λ required for full blink equivalence.

FIG. 13.

Volume and periodicity detection under the cFPLM+ condition with S = 7 × 10−6, V0 = 2.576, m = 4, h0 = 13.

FIG. 13.

Volume and periodicity detection under the cFPLM+ condition with S = 7 × 10−6, V0 = 2.576, m = 4, h0 = 13.

Close modal

We now turn to the dynamics of the OST, which is measured at the surface of the tear film.55,59 Li and Braun16 recovered observed OST cooling rates using an assumed polynomial initial condition that was consistent with observed profiles. We first compute two blink cycles with 5 s interblinks, and upon subsequent opening we compute for a 25 s interblink. The longer interblink is meant to imitate the experimental conditions in which the subjects were asked to keep their eyes open as long as possible.59 Figure 14 shows the computed tear film surface temperature T(x, h, t) at several times starting with 2 s after the eye is completely open, and ending 25 s after. The coolest region is near the middle of the film, which we assume to be the same as the GCC, throughout the open phase, and it approaches the GCC from below as time increases in the model. The location of the minimum temperature inferior to the GCC agrees with experiment.59 This may be explained by the asymmetric eyelid coverage time of the film during a complete blink cycle: the tear film below GCC has more contact time with cool air outside than that above the GCC due to the realistic lid motion during the upstroke and downstroke phase. Efron et al.59 suggested that the asymmetry of the lid margins (the superior lid margin is generally closer to the GCC than the inferior lid margin) contributed to the location of the minimum temperature as well.

FIG. 14.

Surface temperature T(x, h(x, t), t) during the 25 s open phase with S = 7 × 10−6, Bi = 0.0009, he = 1, V0 = 2.576, λ = 0.1, m = 8, and h0 = 25. Two normal blinks (interblink to = 5 s) are applied first to simulate the initial condition for the experiment. The measured blink in the experiment has an interblink of duration to = 25  s.

FIG. 14.

Surface temperature T(x, h(x, t), t) during the 25 s open phase with S = 7 × 10−6, Bi = 0.0009, he = 1, V0 = 2.576, λ = 0.1, m = 8, and h0 = 25. Two normal blinks (interblink to = 5 s) are applied first to simulate the initial condition for the experiment. The measured blink in the experiment has an interblink of duration to = 25  s.

Close modal

In Figure 15 we compute the GCC temperature to compare with the measurement of Efron et al.59 taken 2 s after the eye is completely opened. It is observed that the GCC temperature drops at a greater rate for the first several seconds before flattening out at long times. The simulation predicts the average cooling rate to be 0.045 °C/s. This value lies within the range of the mean cooling rate (0.033 ± 0.024 °C/s) from experiments.59 

FIG. 15.

GCC temperature T(0, h(0, t), t) compared with the experiment data. Here S = 7 × 10−6, he = 1, V0 = 2.576, λ = 0.1, m = 8, and h0 = 25. The initial temperature profile is obtained by simulating two normal blink cycles with 5 s interblink. The figure shows the GCC temperature of the eye blink with 20 s interblink, compared with experiment (from two different subjects in their Figure 3).59 

FIG. 15.

GCC temperature T(0, h(0, t), t) compared with the experiment data. Here S = 7 × 10−6, he = 1, V0 = 2.576, λ = 0.1, m = 8, and h0 = 25. The initial temperature profile is obtained by simulating two normal blink cycles with 5 s interblink. The figure shows the GCC temperature of the eye blink with 20 s interblink, compared with experiment (from two different subjects in their Figure 3).59 

Close modal

We now briefly discuss how the Biot number Bi affects the cooling rate. Bi gives a simple index of the ratio of the heat transfer resistances inside of and at the surface of the tear film. To obtain good agreement with measured cooling rates, Bi = 0.0009 was used in Li and Braun;16 we followed this choice in the results above. The relatively small Bi implies that the heat conduction inside the film is much faster than the heat convection away from its surface. We varied Bi by setting Bi = 0 or increasing Bi by a factor of 5 to Bi = 0.0045.

As can be seen from Figure 16, the GCC cools more quickly with greater Bi. The larger value of Bi can result in physiologically observed temperature drops59 in a 5 s interblink, but those in vivo temperature changes occurs over longer interblinks. This results lends support to our choice of Bi as suitable for matching experimental results in a controlled environment.

FIG. 16.

The surface temperature T(x, h(x, t), t) variation over an interblink (5 s). Curves with different colors correspond to different values of Bi: solid curves (Bi = 0 with GCC cooling rate 0.02 °C/s), dashed curves (Bi = 0.0009 with GCC cooling rate 0.18 °C/s), dashed-dotted curves (Bi = 0.0045 with GCC cooling rate 0.588 °C/s).

FIG. 16.

The surface temperature T(x, h(x, t), t) variation over an interblink (5 s). Curves with different colors correspond to different values of Bi: solid curves (Bi = 0 with GCC cooling rate 0.02 °C/s), dashed curves (Bi = 0.0009 with GCC cooling rate 0.18 °C/s), dashed-dotted curves (Bi = 0.0045 with GCC cooling rate 0.588 °C/s).

Close modal

A normal tear film may have a thinning rate 2.5  μm/min; King-Smith and co-workers73 measured the mean rate of film thinning to be 3.79 ± 4.20  μm/min. In this paper, we have chosen

$\bar{K}$
K¯ to reflect a thinning rate of 3 μm/min. In Figure 17(a), we plot the film thickness at various times during an interblink of 20 s, with two normal blinks in advance to simulate an initial condition. Patients with dry eye syndrome typically show a higher thinning rate than ordinary people. Assuming a thinning rate of 6 μm/min and adjusting
$\bar{K}$
K¯
accordingly, we obtain the result shown in Figure 17(b). Small deviations in the thinning rate from the nominal value occur because of the temperature variations. Both figures suggest that there is a variation in the film thinning rate across the film, agreeing with experimental observations.73 

FIG. 17.

Results for two different thinning rates during a 25 s open phase (matching experiment) with S = 7 × 10−6, he = 1, V0 = 2.576, λ = 0.1, m = 8, and h0 = 25. The open phase shown here was started after two blink cycles with to = 5 s interblinks (open phases) to set up the initial condition for comparison with experimental results.

FIG. 17.

Results for two different thinning rates during a 25 s open phase (matching experiment) with S = 7 × 10−6, he = 1, V0 = 2.576, λ = 0.1, m = 8, and h0 = 25. The open phase shown here was started after two blink cycles with to = 5 s interblinks (open phases) to set up the initial condition for comparison with experimental results.

Close modal

In this paper we simulated tear film dynamics coupled with heat transfer from within the cornea over several blink cycles. We explored this system with both sinusoidal lid motion31 and realistic lid motion,32 using three different flux boundary conditions: no-flux, FPLM and FPLM+. Evaporation was included in the model, so the deformation of the tear film results from a combined effect of capillary thinning and evaporation. All of the computations found that a valley is formed near each ends of the film, corresponding to the phenomenon known as “black line.” This name comes from the experimental observation that thinning is most significant along the line-shaped region near both ends of the film and manifests itself as a dark region when measured with fluorophotometry. With realistic motion and the FPLM+ condition, we computed the thinning rate at all locations across the tear film. We assume 3  μm/min and 6  μm/min average thinning rates to represent the normal and malfunctioning film thinning rates, respectively. Simulated thinning rates agree well with our assumed averages, but the results show that the thinning rate varies across the film: the most rapid thinning is near the black lines due to capillary forces, while away from the black lines the thinning is more uniform.

Full blink equivalence was first suggested and studied by Braun and King-Smith with sinusoidal lid motion.31 Heryudono et al.32 explored this phenomenon further with realistic lid motion and FPLM+ boundary conditions. The numerical results show the necessity of compensating for evaporation with additional influx to obtain the cFPLM (or cFPLM+) condition.24,32 We studied the full blink equivalence with evaporation, where we adopted the compensated FPLM+ condition to restore the volume lost from evaporation. We observed that while λ ≈ 1/6 to 1/8 gives periodic dynamics for sinusoidal motion,31 λ ≈ 0.1 is required for realistic motion to reconstruct periodic behavior, which according to our definition, is effectively a full blink. We conclude that with sufficient flux entering the film to compensate for the evaporation loss, normal blinking eyes could not achieve periodic dynamics without a full closure, while for frequently blinking eyes (sinusoidal lid motion), the eye may not have to be completely closed.

The corneal temperature was explored with both sinusoidal lid motion and realistic lid motion. Craig et al.60 used ocular thermography, a variety of tear physiology tests and measurement of evaporation rates in both control and dry eye subjects in order to correlate physiological factors in dry eye with OST and evaporation from the tear film. They found that there was an increased rate of evaporation in dry eye, as found in a number of previous studies.78–80 They found an evaporation rate of 1.9 × 10−9  g/cm2/s in normals and 4.1 × 10−8  g/cm2/s in dry eye subjects; using their goggle methods showed a factor of 20 increase in some dry eye cases. Still, those goggle methods may give smaller values for the evaporation rate than those measured in other ways.81–83 Dry eye subjects were found to have significantly lower GCC temperature than the controls, and this could be attributed to increased evaporation. Also, the rate of cooling in the cornea was significantly higher in dry eyes, which may be caused by increased evaporation as well. Notably, Craig et al.60 did not find any correlation between reduced non-invasive break up time (time to first rupture) and evaporation. Two studies have found that evaporation can be reduced in dry eye,84,85 but these authors point out that their subjects are aqueous deficient dry eyes, and Craig et al.60 speculate that this may be responsible for the reduced evaporation rate.

Both temperature profiles (Figures 7 and 11) correlate directly with the lid motion, where a cool island is formed at the end of every blink cycle and serves as the initial condition for the next blink. We can also observe that a nearly periodic state for the cornea temperature is achieved within only one or two blink cycles. Since for the realistic lid motion the eye is kept open for a much longer period relative to the blink than in the sinusoidal case, the temperature decrease extends well into the cornea. In experiments56,59 ocular surface temperature is measured using infrared technology, and it is observed that the coolest point lies slightly inferior to the GCC. Our computational results captured this phenomenon using the simulated initial condition from two successive normal blinks (5 s interblinks). The transient cooling of the eye appears to contribute to the minimum temperature being inferior to the GCC, in addition to the unbalanced lid coverage during a blink cycle. GCC temperature computed based on this initial condition agrees well with experiment,59 and the cooling rate 0.045 °C/s computed lies in the range of mean cooling rate (0.033 ± 0.024 °C/s).

Heat transfer from the posterior side of the tear film was incorporated into a model for tear film dynamics. Evaporation of water from the tear film and a Newton law of cooling approximating convection outside the tear film are the mechanisms for heat loss from the anterior surface of the tear film into the air. In the model, diffusion of heat within the eye and the tear film is the only mechanism of heat transfer in those areas. The dynamics of the temperature field shows periodicity, with a small cool region trapped inside the eye during a blink. The coolest location during the interblink is just inferior to the GCC, and the minimum is within the experimentally observed range. Our results suggest that the transient thermal effects due to blinking help locate the minimum temperature below the GCC, not simply the open lid margin locations alone as suggested by Efron et al.59 

The tear film dynamics with regard to periodic eyelid motion are similar to the dynamics computed in Heryudono et al.32 where evaporation is absent. In that work, for realistic lid motion, the solutions became periodic when λ ⩽ 0.13; in the current model, with evaporation, a compensating influx of fluid, and realistic lid motion, we find that only λ ≈ 0.1 (a full blink) results in periodic solutions. However, some care is needed to compensate for the evaporative mass loss in order seek periodic solutions in the first place.

Future directions include incorporating two-dimensional tear films on eye-opening-shaped domains that approximate the tear film to be solved together with three dimensional regions of heat diffusion that approximate the underlying eye. Such models may use different properties for cornea/aqueous humor substrates as in this paper than for conjunctival/stroma/vitreous substrates that are posterior to the tear film outside the corneal region. Furthermore, more complex treatments of the tear film, e.g., with separate lipid and aqueous layers with evaporation,36,86 may add insights into tear film and OST dynamics, and the OST may be an important input into the properties of the lipid layer.45 

This material is based on work supported by the National Science Foundation under Grant No. DMS-1022706. The authors thank Mr. Longei Li for helpful conversations and input data.

The end motion function that we refer to as realistic lid motion is as follows:

\begin{equation}X(t) = \left\lbrace \begin{array}{l l l}1-2\lambda -2(1-\lambda )\left(\frac{t}{\Delta t_{\text{co}}}\right)^2\exp {\left[1-\left(\frac{t}{\Delta t_{\text{co}}}\right)^2\right]}, & \ 0 \le t\le \Delta t_{\text{co}} \\-1, & \ \Delta t_{\text{co}} \le t \le \Delta t_{\text{co}}+\Delta t_{\text{o}}\\-1 + 2(1-\lambda )\left(\frac{t-\Delta t_{\text{co}} - \Delta t_{\text{o}}}{\Delta t_{\text{oc}}}\right)^2\exp {\left[1-\left(\frac{t-\Delta t_{\text{co}} - \Delta t_{\text{o}}}{\Delta t_{\text{oc}}}\right)^2\right]}. & \ \Delta t_{\text{co}}+\Delta t_{\text{o}} \le t\le \Delta t_{\text{bc}} \\\end{array} \right.\end{equation}
X(t)=12λ2(1λ)tΔtco2exp1tΔtco2,0tΔtco1,ΔtcotΔtco+Δto1+2(1λ)tΔtcoΔtoΔtoc2exp1tΔtcoΔtoΔtoc2.Δtco+ΔtotΔtbc
(A1)

Here Δtbc = Δtco + Δto + Δtoc is the non-dimensionalized period of a complete normal blink cycle.

Besides the fluxes that are caused solely by lid motion, Heryudono et al.32 modified the FPLM conditions by adding more terms representing influx from the lacrimal gland (subscript lg) and outflux via the puncta (subscript p), respectively. The flux functions are then as follows:

\begin{eqnarray}Q_{\text{top}} = -X_t h_e/2 -2 f_{\text{out}} Q_{\text{0p}} \exp {\left[-\left(\frac{t-t_{\text{out}}}{\Delta t_{\text{p}}}\right)^2\right]} \nonumber \\+\, f_{\text{top}} Q_{\text{0lg}}\exp {\left[-\left(\frac{t-t_{\text{in}}}{\Delta t_{\text{co}}/2}\right)^2\right]},\qquad\qquad\end{eqnarray}
Qtop=Xthe/22foutQ0pexpttoutΔtp2+ftopQ0lgexpttinΔtco/22,
(A2)
\begin{eqnarray}Q_{\text{bot}} = -&2 (1-f_{\text{out}}) Q_{\text{0p}} \exp {\left[-\left(\frac{t-t_{\text{out}}}{\Delta t_{\text{p}}}\right)^2\right]} \nonumber \\+& (1-f_{\text{top}}) Q_{\text{0lg}}\exp {\left[-\left(\frac{t-t_{\text{in}}}{\Delta t_{\text{co}}/2}\right)^2\right]}.\end{eqnarray}
Qbot=2(1fout)Q0pexpttoutΔtp2+(1ftop)Q0lgexpttinΔtco/22.
(A3)

Here Q0lg represents the supply of tear fluid from the lacrimal gland and Q0p represents the drainage through the puncta. We refer to the flux conditions that use these functions as the FPLM+ condition, and use the same parameters as specified in Heryudono et al.32 unless otherwise noted. (See Table III.)

Table III.

Parameter values in the definitions of lid motion and FPLM+ flux condition.32 

Nondimensionalized parameters for lid motion and flux condition
ParameterDescriptionValue
Δtco Duration of upstroke 0.176 
Δto Duration of interblink 
Δtoc Duration of downstroke 0.88 
Δtp Duration of punctal drainage 0.5 
Δtlg Duration of lacrimal gland supply 0.088 
tin Time location of influx peak 0.176 
tout Time location of outflux peak 1.08 
ftop Fraction of supply/drainage from the superior side 0.7 
fbot Fraction of supply/drainage from the inferior side 0.3 
Qlg Tear supply from lacrimal gland 0.337 
Qp Drainage due to the puncta 0.0297 
Nondimensionalized parameters for lid motion and flux condition
ParameterDescriptionValue
Δtco Duration of upstroke 0.176 
Δto Duration of interblink 
Δtoc Duration of downstroke 0.88 
Δtp Duration of punctal drainage 0.5 
Δtlg Duration of lacrimal gland supply 0.088 
tin Time location of influx peak 0.176 
tout Time location of outflux peak 1.08 
ftop Fraction of supply/drainage from the superior side 0.7 
fbot Fraction of supply/drainage from the inferior side 0.3 
Qlg Tear supply from lacrimal gland 0.337 
Qp Drainage due to the puncta 0.0297 

We define the pressure p(x, t) = −(Shxx + Ah−3) as a new dependent variable. The thin film equations (14) and (15) then become

\begin{equation}q = -\frac{h^3}{12}p_x , \ \ J = \left[\frac{T_0 + \text{Bi} T_{\infty } h}{1+ \text{Bi} h} - \delta p\right] \left[\bar{K} + \frac{h}{1+\text{Bi} h}\right]^{-1}.\end{equation}
q=h312px,J=T0+BiTh1+BihδpK¯+h1+Bih1.
(B1)

Along with the conservation statement (13), these four equations define a closed system.

Since spectral collocation methods lead to dense matrices and the computational cost to solve them grows as the third power of the number of unknowns, we used Chebyshev meshes with different density for X(t) ⩽ x ⩽ 1, y = 0 (the 1D film domain) and

$-1 \le x \le 1, 0 \le \tilde{y} \le \tilde{y}_c$
1x1,0ỹỹc (the 2D temperature domain posterior to the film). We used a fine grid for the film domain (200 Chebyshev points) and a coarse grid for the temperature domain (20 × 20 Chebyshev points in the x and y directions). At each time step, we used barycentric interpolation87 to communicate between the two meshes at the cornea/tear film interface X(t) ⩽ x ⩽ 1 along y = 0.

The film domain X(t) ⩽ x ⩽ 1 is mapped to a fixed domain −1 ⩽ ξ ⩽ 1 via the time-varying coordinate change

\begin{equation}\xi = 1 - 2\frac{1-x}{1-X(t)}.\end{equation}
ξ=121x1X(t).
(B2)

Using the change of variables h(x, t) = H(ξ(t), t) and q(x, t) = Q(ξ(t), t), we rewrite Eqs. (13)–(29) as

\begin{equation}H_t = \frac{1-\xi }{1-X} X_t H_{\xi } - EJ -\left(\frac{2}{1-X}\right)Q_{\xi },\end{equation}
Ht=1ξ1XXtHξEJ21XQξ,
(B3)
\begin{eqnarray}\hspace*{5pt}Q = &\frac{H^3}{12} \left(\frac{2}{1-X}\right) \left[S \left(\frac{2}{1-X}\right)^2 H_{\xi \xi } + A H^{-3}\right]_{\xi } \nonumber \\&+ X_t \frac{1-\xi }{2}\frac{H}{2},\end{eqnarray}
Q=H31221XS21X2Hξξ+AH3ξ+Xt1ξ2H2,
(B4)
\begin{eqnarray}J = &\left[ \frac{T_0 + \text{Bi} T_\infty H}{1+ \text{Bi} H} - \delta \left( S \left(\frac{2}{1-X}\right)^2 H_{\xi \xi } + A_1 H^{-3} \right) \right] \left[ K + \frac{H}{1 + \text{Bi} H} \right]^{-1}.\end{eqnarray}
J=T0+BiTH1+BiHδS21X2Hξξ+A1H3K+H1+BiH1.
(B5)

The initial condition becomes

\begin{equation}H(\xi ,0) = H_{\text{min}} +(h_0 - H_{\text{min}}) \xi ^{m},\end{equation}
H(ξ,0)=Hmin+(h0Hmin)ξm,
(B6)

and the boundary conditions become

\begin{eqnarray}H(\pm 1,t) = h_0, \ \ Q(1,t) = -Q_{\text{bot}} , \nonumber \\Q(-1,t) = X_t h_0 + Q_{\text{top}}.\end{eqnarray}
H(±1,t)=h0,Q(1,t)=Qbot,Q(1,t)=Xth0+Qtop.
(B7)

Once the evolution equations are discretized in space, the result is a finite-dimensional system of DAEs. We used ode15s in MATLAB (similar to the backward differentiation formula) to solve the discrete system.

1.
F. J.
Holly
and
M. A.
Lemp
, “
Tear physiology and dry eyes
,”
Rev. Surv. Ophthalmol.
22
,
69
(
1977
).
2.
S.
Mishima
, “
Some physiological aspects of the precorneal tear film
,”
Arch. Ophthalmol.
73
,
233
(
1965
).
3.
N.
Ehlers
, “
The precorneal film: Biomicroscopical, histological and chemical investigations
,”
Acta Ophthalmol. Suppl.
81
,
3
(
1965
).
4.
M. S.
Norn
, “
Semiquantitative interference study of fatty layer of precorneal film
,”
Acta Ophthalmol.
57
,
766
(
1979
).
5.
A. J.
Bron
,
J. M.
Tiffany
,
S. M.
Gouveia
,
N.
Yokoi
, and
L. W.
Voon
, “
Functional aspects of the tear film lipid layer
,”
Exp. Eye Res.
78
,
347
(
2004
).
6.
P. E.
King-Smith
,
B. A.
Fink
,
J. J.
Nichols
,
K. K.
Nichols
,
R. J.
Braun
, and
G. B.
McFadden
, “
The contribution of lipid layer movement to tear film thinning and breakup
,”
Invest. Ophthalmol. Visual Sci.
50
,
2747
(
2009
).
7.
H.-B.
Chen
,
S.
Yamabayashi
,
B.
Ou
,
Y.
Tanaka
, and
S.
Ohno
, “
Structure and composition of rat precorneal tear film: A study by in vivo cryofixation
,”
Invest. Ophthalmol. Visual Sci.
38
,
381
(
1997
).
8.
I. K.
Gipson
, “
Distribution of mucins at the ocular surface
,”
Exp. Eye Res.
78
,
379
(
2004
).
9.
B.
Govindarajan
and
I. K.
Gipson
, “
Membrane-tethered mucins have multiple functions on the ocular surface
,”
Exp. Eye Res.
90
,
655
(
2010
).
10.
P. E.
King-Smith
,
B. A.
Fink
,
R. M.
Hill
,
K. W.
Koelling
, and
J. M.
Tiffany
, “
The thickness of the tear film
,”
Curr. Eye Res.
29
,
357
(
2004
).
11.
P. E.
King-Smith
,
B. A.
Fink
,
J. J.
Nichols
,
K. K.
Nichols
, and
R. M.
Hill
, “
Interferometric imaging of the full thickness of the precorneal tear film
,”
J. Opt. Soc. Am. A
23
,
2097
(
2006
).
12.
J.
Wang
,
D.
Fonn
,
T. L.
Simpson
, and
L.
Jones
, “
Precorneal and pre- and postlens tear film thickness measured indirectly with optical coherence tomography
,”
Invest. Ophthalmol. Visual Sci.
44
,
2524
(
2003
).
13.
J. R.
Palakuru
,
J.
Wang
, and
J. V.
Aquavella
, “
Effect of blinking on tear dynamics
,”
Invest. Ophthalmol. Visual Sci.
48
,
3032
(
2007
).
14.
M. E.
Johnson
and
P. J.
Murphy
, “
Temporal changes in the tear menisci following a blink
,”
Exp. Eye Res.
83
,
517
(
2006
).
15.
W. W.
Harrison
,
C. G.
Begley
,
H.
Lui
,
M.
Chen
,
M.
Garcia
, and
J. A.
Smith
, “
Menisci and fullness of the blink in dry eye
,”
Optom. Visual Sci.
85
,
706
(
2008
).
16.
L.
Li
and
R. J.
Braun
, “
A model for the human tear film with heating from within the eye
,”
Phys. Fluids.
24
,
062103
(
2012
).
17.
R. J.
Braun
, “
Dynamics of the tear film
,”
Annu. Rev. Fluid Mech.
44
,
267
(
2012
).
18.
R. J.
Braun
,
R.
Usha
,
G. B.
McFadden
,
T. A.
Driscoll
,
L. P.
Cook
, and
P. E.
King-Smith
, “
Thin film dynamics on a prolate spheroid with application to the cornea
,”
J. Eng. Math.
73
,
121
(
2012
).
19.
K. L.
Miller
,
K. A.
Polse
, and
C. J.
Radke
, “
Black line formation and the “perched” human tear film
,”
Curr. Eye Res.
25
,
155
(
2002
).
20.
H.
Wong
,
I.
Fatt
, and
C. J.
Radke
, “
Deposition and thinning of the human tear film
,”
J. Colloid Interface Sci.
184
,
44
(
1996
).
21.
A.
Sharma
,
S.
Tiwari
,
R.
Khanna
, and
J. M.
Tiffany
, “
Hydrodynamics of meniscus-induced thinning of the tear film
,” in
Lacrimal Gland, Tear Film, and Dry Eye Syndromes 2
,
Advances in Experimental Medicine and Biology
Vol.
438
, edited by
D. A.
Sullivan
,
D. A.
Dartt
, and
M. A.
Meneray
(
Springer
,
1998
), p.
425
.
22.
J. E.
McDonald
and
S.
Brubaker
, “
Meniscus-induced thinning of tear films
,”
Am. J. Ophthalmol.
72
,
139
(
1971
).
23.
R. E.
Berger
and
S.
Corrsin
, “
A surface tension gradient mechanism for driving the pre-corneal tear film after a blink
,”
J. Biomech.
7
,
225
(
1974
).
24.
M. B.
Jones
,
C. P.
Please
,
D. L. S.
McElwain
,
G. R.
Fulford
,
A. P.
Roberts
, and
M. J.
Collins
, “
Dynamics of tear film deposition and drainage
,”
Math. Med. Biol.
22
,
265
(
2005
).
25.
M. B.
Jones
,
D. L. S.
McElwain
,
G. R.
Fulford
,
M. J.
Collins
, and
A. P.
Roberts
, “
The effect of the lipid layer on tear film behavior
,”
Bull. Math. Biol.
68
,
1355
(
2006
).
26.
E.
Aydemir
,
C. J. W.
Breward
, and
T. P.
Witelski
, “
The effect of polar lipids on tear film dynamics
,”
Bull. Math. Biol.
73
,
1171
(
2010
).
27.
K. N.
Winter
,
D. M.
Anderson
, and
R. J.
Braun
, “
A model for wetting and evaporation of a post-blink precorneal tear film
,”
Math. Med. Biol.
27
,
211
(
2010
).
28.
R. J.
Braun
and
A. D.
Fitt
, “
Modeling the drainage of the precorneal tear film after a blink
,”
Math. Med. Biol.
20
,
1
(
2003
).
29.
Y. L.
Zhang
,
O. K.
Matar
, and
R. V.
Craster
, “
Analysis of tear film rupture: Effect of non-Newtonian rheology
,”
J. Colloid Interface Sci.
262
,
130
(
2003
).
30.
Y. L.
Zhang
,
O. K.
Matar
, and
R. V.
Craster
, “
Rupture analysis of the corneal mucus layer of the tear film
,”
Mol. Simul.
30
,
167
(
2004
).
31.
R. J.
Braun
and
P. E.
King-Smith
, “
Model problems for the tear film in a blink cycle: Single equation models
,”
J. Fluid Mech.
586
,
465
(
2007
).
32.
A.
Heryudono
,
R. J.
Braun
,
T. A.
Driscoll
,
L. P.
Cook
,
K. L.
Maki
, and
P. E.
King-Smith
, “
Single-equation models for the tear film in a blink cycle: Realistic lid motion
,”
Math. Med. Biol.
24
,
347
(
2007
).
33.
K. L.
Maki
,
R. J.
Braun
,
T. A.
Driscoll
, and
P. E.
King-Smith
, “
An overset grid method for the study of reflex tearing
,”
Math. Med. Biol.
25
,
187
(
2008
).
34.
L.
Jossic
,
P.
Lefevre
,
C.
de Loubens
,
A.
Magnin
, and
C.
Corre
, “
The fluid mechanics of shear-thinning tear substitutes
,”
J. Non-Newton. Fluid Mech.
161
,
1
(
2009
).
35.
V.
Zubkov
,
C. J. W.
Breward
, and
E. A.
Gaffney
, “
Coupling fluid and solute dynamics within the ocular surface tear film: A modelling study of black line osmolarity
,”
Bull. Math. Biol.
74
,
2062
(
2012
).
36.
M.
Bruna
and
C. J. W.
Breward
, “
The influence of nonpolar lipids on tear film dynamics
,”
J. Fluid Mech.
746
,
565
605
(
2014
).
37.
K. L.
Maki
,
R. J.
Braun
,
P.
Ucciferro
,
W. D.
Henshaw
, and
P. E.
King-Smith
, “
Tear film dynamics on an eye-shaped domain II: Flux boundary conditions
,”
J. Fluid Mech.
647
,
361
(
2010
).
38.
K. L.
Maki
,
R. J.
Braun
,
W. D.
Henshaw
, and
P. E.
King-Smith
, “
Tear film dynamics on an eye-shaped domain I: Pressure boundary conditions
,”
Math. Med. Biol.
27
,
227
(
2010
).
39.
M. G.
Doane
, “
Interactions of eyelids and tears in corneal wetting and the dynamics of the normal human eyeblink
,”
Am. J. Ophthalmol.
89
,
507
(
1980
).
40.
M.
Lorber
, “
Gross characteristics of normal human lacrimal glands
,”
Ocul. Surf.
5
,
13
(
2007
).
41.
D. M.
Maurice
, “
The dynamics and drainage of tears
,”
Int. Ophthalmol. Clin.
13
,
103
(
1973
).
42.
J. M.
Tiffany
, “
The viscosity of human tears
,”
Int. Ophthalmol.
15
,
371
(
1991
).
43.
J. C.
Pandit
,
B.
Nagyová
,
A. J.
Bron
, and
J. M.
Tiffany
, “
Physical properties of stimulated and unstimulated tears
,”
Exp. Eye Res.
68
,
247
(
1999
).
44.
D. L.
Leiske
,
S. R.
Raju
,
H. A.
Ketelson
,
T. J.
Millar
, and
G. G.
Fuller
, “
The interfacial viscoelastic properties and structures of human and animal Meibomian lipids
,”
Exp. Eye Res.
90
,
598
(
2010
).
45.
D. L.
Leiske
,
C. I.
Leiske
,
D. R.
Leiske
,
M. F.
Toney
,
M.
Senchyna
,
H. A.
Ketelson
,
D. L.
Meadows
, and
G. G.
Fuller
, “
Temperature-induced transitions in the structure and interfacial rheology of human meibum
,”
Biophys. J.
102
,
369
(
2012
).
46.
J. P.
McCulley
and
W.
Shine
, “
A compositional based model for the tear film lipid layer
,”
Trans. Am. Ophthalmol. Soc.
95
,
79
(
1997
).
47.
B.
Nagyová
and
J. M.
Tiffany
, “
Components of tears responsible for surface tension
,”
Curr. Eye Res.
19
,
4
(
1999
).
48.
H.
Owens
and
J.
Phillips
, “
Spread of the tears after a blink: Velocity and stabilization time in healthy eyes
,”
Cornea
20
,
484
(
2001
).
49.
P.
Mudgil
,
M.
Torres
, and
T. J.
Millar
, “
Adsorption of lysozyme to phospholipid and meibomian lipid monolayer films
,”
Colloids Surf., B
48
,
128
(
2006
).
50.
P.
Mudgil
and
T. J.
Millar
, “
Adsorption of apo- and holo-tear lipocalin to a bovine Meibomian lipid film
,”
Exp. Eye Res.
86
,
622
(
2008
).
51.
V. S.
Ajaev
, “
Spreading of thin volatile liquid droplets on uniformly heated surfaces
,”
J. Fluid Mech.
528
,
279
(
2005
).
52.
V. S.
Ajaev
, “
Evolution of dry patches in evaporating liquid films
,”
Phys. Rev. E
72
,
031605
(
2005
).
53.
V. S.
Ajaev
and
G. M.
Homsy
, “
Steady vapor bubbles in rectangular microchannels
,”
J. Colloid Interface Sci.
240
,
259
(
2001
).
54.
C.
Sodtke
,
V. S.
Ajaev
, and
P.
Stephan
, “
Dynamics of volatile liquid droplets on heated surfaces: theory versus experiment
,”
J. Fluid Mech.
610
,
343
(
2008
).
55.
C.
Purslow
and
J. S.
Wolffsohn
, “
Ocular surface temperature: A review
,”
Eye Contact Lens
31
,
117
(
2005
).
56.
R.
Mapstone
, “
Measurement of corneal temperature
,”
Exp. Eye Res.
7
,
237
(
1968
).
57.
R.
Mapstone
, “
Determinants of corneal temperature
,”
Br. J. Ophthalmol.
52
,
729
(
1968
).
58.
R.
Mapstone
, “
Normal thermal patterns in cornea and periorbital skin
,”
Br. J. Ophthalmol.
52
,
818
(
1968
).
59.
N.
Efron
,
G.
Young
, and
N. A.
Brennan
, “
Ocular surface temperature
,”
Curr. Eye Res.
8
,
901
(
1989
).
60.
J. P.
Craig
,
I.
Singh
,
A.
Tomlinson
,
P. B.
Morgan
, and
N.
Efron
, “
The role of tear physiology in ocular surface temperature
,”
Eye
14
,
635
(
2000
).
61.
C.
Purslow
,
J. S.
Wolffsohn
, and
J.
Santodomingo-Rubido
, “
The effect of contact lens wear on dynamic ocular surface temperature
,”
Contact Lens Anterior Eye
28
,
29
(
2005
).
62.
J.
Aliò
and
M.
Padron
, “
Influence of age on the temperature of the anterior segment of the eye
,”
Ophthalmic Res.
14
,
153
(
1982
).
63.
J. A.
Scott
, “
A finite element model of heat transport in the human eye
,”
Phys. Med. Biol.
33
,
227
(
1988
).
64.
J. A.
Scott
, “
The computation of temperature rises in the human eye induced by infrared radiation
,”
Phys. Med. Biol.
33
,
243
(
1988
).
65.
E. Y. K.
Ng
and
E. H.
Ooi
, “
FEM simulation of the eye structure with bioheat analysis
,”
Comput. Methods Prog. Biomed.
82
,
268
(
2006
).
66.
E. Y. K.
Ng
and
E. H.
Ooi
, “
Ocular surface temperature: A 3D FEM prediction using bioheat equation
,”
Comput. Biol. Med.
37
,
829
(
2007
).
67.
J. J.
Heys
and
V. H.
Barocas
, “
A Boussinesq model of natural convection in the human eye and the formation of Krukenberg's spindle
,”
Ann. Biomed. Eng.
30
,
392
(
2002
).
68.
E. H.
Ooi
and
E. Y. K.
Ng
, “
Simulation of aqueous humor hydrodynamics in human eye heat transfer
,”
Comput. Biol. Med.
38
,
252
(
2008
).
69.
A.
Karampatzakis
and
T.
Samaras
, “
Numerical model of heat transfer in the human eye with consideration of fluid dynamics of the aqueous humour
,”
Phys. Med. Biol.
55
,
5653
(
2010
).
70.
C. R.
Canning
,
M. J.
Greaney
,
J. N.
Dewynne
, and
A. D.
Fitt
, “
Fluid flow in the anterior chamber of the eye
,”
IMA J. Maths. Appl. Med. Biol.
19
,
31
(
2002
).
71.
A. D.
Fitt
and
G.
Gonzalez
, “
Fluid mechanics of the human eye: Aqueous humour flow in the anterior chamber
,”
Bull. Math. Biol.
68
,
53
(
2006
).
72.
J. H.
Siggers
and
C. R.
Ethier
, “
Fluid Mechanics of the Eye
,”
Annu. Rev. Fluid Mech.
44
,
347
(
2012
).
73.
J. J.
Nichols
,
G. L.
Mitchell
, and
P. E.
King-Smith
, “
Thinning rate of the precorneal and prelens tear films
,”
Invest. Ophthalmol. Visual Sci.
46
,
2353
(
2005
).
74.
A.
Oron
,
S. H.
Davis
, and
S. G.
Bankoff
, “
Long-scale evolution of thin liquid films
,”
Rev. Mod. Phys.
69
,
931
(
1997
).
75.
R. V.
Craster
and
O. K.
Matar
, “
Dynamics and stability of thin liquid films
,”
Rev. Mod. Phys.
81
,
1131
(
2009
).
76.
A.
Berke
and
S.
Mueller
, “
Einfluss des lidschlages auf die Kontaktlinse und die zugrundeliegenden Kräfte
,”
die Kontaktlinse
1
,
17
(
1996
).
77.
L. N.
Trefethen
,
Spectral Methods in Matlab
(
SIAM
,
Philadelphia
,
2000
).
78.
M.
Rolando
and
M. F.
Refojo
, “
Tear evaporimeter for measuring water evaporation rate from the tear film under controlled conditions in humans
,”
Exp. Eye Res.
36
,
25
(
1983
).
79.
W. D.
Mathers
,
G.
Binarao
, and
M.
Petroll
, “
Ocular water evaporation and the dry eye: a new measuring device
,”
Cornea
12
,
335
(
1993
).
80.
W. D.
Mathers
, “
Ocular evaporation in meibomian gland dysfunction and dry eye
,”
Ophthalmology
100
,
347
(
1993
).
81.
A.
Tomlinson
,
M. G.
Doane
, and
A.
McFayden
, “
Inputs and outputs of the lacrimal system: Review of production and evaporative loss
,”
Ocul. Surf.
7
,
186
(
2009
).
82.
S. H.
Kimball
,
P. E.
King-Smith
, and
J. J.
Nichols
, “
Evidence for the major contribution of evaporation to tear film thinning between blinks
,”
Invest. Opthalmol. Visual Sci.
51
,
6294
(
2010
).
83.
W. D.
Mathers
, “
Evaporation from the ocular surface
,”
Exp. Eye Res.
78
,
389
(
2004
).
84.
H.
Fujishima
,
I.
Toda
,
M.
Yamada
,
N.
Sato
, and
K.
Tsubota
, “
Corneal Temperature in patients with dry eye evaluated by infrared radiation thermometry
,”
Br. J. Ophthalmol.
80
,
29
(
1996
).
85.
A.
Mori
,
Y.
Oguchi
,
Y.
Okusawa
,
M.
Onon
,
H.
Fujishima
, and
K.
Tsubota
, “
Use of high-speed, high-resolution thermography to evaluate the tear film layer
,”
Am. J. Ophthalmol.
124
,
729
(
1997
).
86.
C.-C.
Peng
,
C.
Cerretani
,
R. J.
Braun
, and
C. J.
Radke
, “
Evaporation-driven instability of the precorneal tear film
,”
Adv. Coll. Interface Sci.
206
,
250
264
(
2014
).
87.
J.-P.
Berrut
and
L. N.
Trefethen
, “
Barycentric Lagrange interpolation
,”
SIAM Rev.
46
,
501
(
2004
).