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.

## I. INTRODUCTION

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 film^{2,3} has an anterior layer of lipids that is in contact with air^{4–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 layer^{1} 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 blink^{10–12} and has a considerably thicker meniscus (about 0.1 mm or more) around the lid margins^{13–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 Braun^{16} 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, Doane^{39} 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-dimensional^{32,33} and two-dimensional^{37,38} models.

We briefly comment on the tear film fluid properties. A sample of tears combining all of its components is mildly shear thinning^{42} but weakly elastic.^{43} Recent measurements of meibomian lipids (meibum) alone from various mammals show that meibum may have significant elasticity^{44} 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 lipids^{46,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 Braun^{16} 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

^{16}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 PCTF^{59,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 study^{62} 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.

Scott^{63,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 developed^{63} a detailed axisymmetric FE model of heat transfer within the eye and to the external environment beyond the anterior of the eye. Scott^{64} 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 eye^{65,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 two^{67,68} and three^{69} 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 Ethier^{72} 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 Braun^{16} 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 model^{16} 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.

## II. FORMULATION

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.

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*^{′}) ⩽ *x* ⩽ *L*^{′} 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

*d*

^{′}represents the typical thickness of the tear film, making the substrate thick compared with the tear film.

^{16,54}In physical terms,

*c*

_{p}and thermal conductivity

*k*. Dimensional parameter values are given in Table I.

Dimensional parameters . | ||
---|---|---|

Parameter . | Description . | Value . |

d^{′} | Characteristic thickness of tear film | 5 × 10^{−6} m |

L^{′} | Half length of the palpebral fissure | 5 × 10^{−3} m |

$h^{\prime }_e$ $he\u2032$ | Thickness of tear fluid under eyelid | 0.25d^{′} to d^{′} |

$L^{\prime }_c$ $Lc\u2032$ | Corneal thickness | 5 × 10^{−4} m |

$y^{\prime }_c$ $yc\u2032$ | Model depth posterior to film | $5L^{\prime }_c$ $5Lc\u2032$ |

μ | Viscosity^{42} | 1.3 × 10^{−3} Pa s |

σ | Surface tension^{47} | 0.045 N m^{−1} |

k | Tear film thermal conductivity^{63} (water) | 0.68 W m^{−1} K^{−1} |

c_{p} | Specific heat (water) | 4184 J kg^{−1} (C°)^{−1} |

k_{c} | Corneal thermal conductivity^{63} | 0.58 W m^{−1} K^{−1} |

ρ | Density (water) | 10^{3} kg m^{−3} |

${\cal L}_m$ $Lm$ | Latent heat of vaporization (water) | 2.3 × 10^{6} J kg^{−1} |

$T^{\prime }_s$ $Ts\u2032$ | Saturation temperature (assumed) | 27 °C |

$T^{\prime }_{\text{B}}$ $TB\u2032$ | Body temperature | 37 °C |

g | Gravitational acceleration | 9.81 m s^{−1} |

κ_{c} | Thermal diffusivity of cornea^{63} | 1.9 × 10^{−7} m^{2} s^{−1} |

A* | Hamaker constant^{27} | 3.5 × 10^{−19} Pa m^{3} |

U^{′} | Representative interblink speed | $U^{\prime } = 5\ \textrm {mm\,s}^{-1}$ $U\u2032=5 mm s\u22121$ |

α | Pressure coefficient for evaporation^{27} | 3.6 × 10^{−2} K Pa^{−1} |

K | Non-equilibrium coefficient^{16} | 1.5 × 10^{5} K m^{2} s kg^{−1} |

Dimensional parameters . | ||
---|---|---|

Parameter . | Description . | Value . |

d^{′} | Characteristic thickness of tear film | 5 × 10^{−6} m |

L^{′} | Half length of the palpebral fissure | 5 × 10^{−3} m |

$h^{\prime }_e$ $he\u2032$ | Thickness of tear fluid under eyelid | 0.25d^{′} to d^{′} |

$L^{\prime }_c$ $Lc\u2032$ | Corneal thickness | 5 × 10^{−4} m |

$y^{\prime }_c$ $yc\u2032$ | Model depth posterior to film | $5L^{\prime }_c$ $5Lc\u2032$ |

μ | Viscosity^{42} | 1.3 × 10^{−3} Pa s |

σ | Surface tension^{47} | 0.045 N m^{−1} |

k | Tear film thermal conductivity^{63} (water) | 0.68 W m^{−1} K^{−1} |

c_{p} | Specific heat (water) | 4184 J kg^{−1} (C°)^{−1} |

k_{c} | Corneal thermal conductivity^{63} | 0.58 W m^{−1} K^{−1} |

ρ | Density (water) | 10^{3} kg m^{−3} |

${\cal L}_m$ $Lm$ | Latent heat of vaporization (water) | 2.3 × 10^{6} J kg^{−1} |

$T^{\prime }_s$ $Ts\u2032$ | Saturation temperature (assumed) | 27 °C |

$T^{\prime }_{\text{B}}$ $TB\u2032$ | Body temperature | 37 °C |

g | Gravitational acceleration | 9.81 m s^{−1} |

κ_{c} | Thermal diffusivity of cornea^{63} | 1.9 × 10^{−7} m^{2} s^{−1} |

A* | Hamaker constant^{27} | 3.5 × 10^{−19} Pa m^{3} |

U^{′} | Representative interblink speed | $U^{\prime } = 5\ \textrm {mm\,s}^{-1}$ $U\u2032=5 mm s\u22121$ |

α | Pressure coefficient for evaporation^{27} | 3.6 × 10^{−2} K Pa^{−1} |

K | Non-equilibrium coefficient^{16} | 1.5 × 10^{5} K m^{2} 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

*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

*T*is made non-dimensionalized as

*k*is the thermal conductivity, and

. | Dimensionless parameters . | . |
---|---|---|

ε | $\frac{d^{\prime }}{L^{\prime }}$ $d\u2032L\u2032$ | 10^{−3} |

Re | $\frac{\rho U^{\prime } L^{\prime }}{\mu }$ $\rho U\u2032L\u2032\mu $ | 20 |

Pr | $\frac{\mu c_P}{k}$ $\mu cPk$ | 8.0 |

E | $\frac{k(T^{\prime }_{\text{B}}-T^{\prime }_s)}{{d^{\prime } {\cal L}_m} \epsilon \rho U^{\prime }}$ $k(TB\u2032\u2212Ts\u2032)d\u2032Lm\epsilon \rho U\u2032$ | 118.3 |

S | $\frac{\sigma \epsilon ^3}{\mu U^{\prime }}$ $\sigma \epsilon 3\mu U\u2032$ | 6.92 × 10^{−6} |

$\bar{K}$ $K\xaf$ | $\frac{kK}{{d^{\prime } {\cal L}_m}}$ $kKd\u2032Lm$ | 8.9 × 10^{3} |

h_{e} | $\frac{h^{\prime }_e}{d^{\prime }}$ $he\u2032d\u2032$ | 1 |

δ | $\frac{\alpha \mu U^{\prime }}{\epsilon ^2L(T^{\prime }_{\text{B}}-T^{\prime }_s)}$ $\alpha \mu U\u2032\epsilon 2L(TB\u2032\u2212Ts\u2032)$ | 4.66 |

A | $\frac{A^*}{{L^{\prime } d^{\prime }} \mu U^{\prime }}$ $A*L\u2032d\u2032\mu U\u2032$ | 2.14 × 10^{−6} |

$\tilde{k}$ $k\u0303$ | k/k_{c} | 1.17 |

P_{T} | $\frac{L^{\prime } \kappa _c}{U^{\prime } {{L^{\prime }_c}^2}}$ $L\u2032\kappa cU\u2032Lc\u20322$ | 0.76 |

Bi | $\frac{h_cd^{\prime }}{k}$ $hcd\u2032k$ | 0-0.1 |

. | Dimensionless parameters . | . |
---|---|---|

ε | $\frac{d^{\prime }}{L^{\prime }}$ $d\u2032L\u2032$ | 10^{−3} |

Re | $\frac{\rho U^{\prime } L^{\prime }}{\mu }$ $\rho U\u2032L\u2032\mu $ | 20 |

Pr | $\frac{\mu c_P}{k}$ $\mu cPk$ | 8.0 |

E | $\frac{k(T^{\prime }_{\text{B}}-T^{\prime }_s)}{{d^{\prime } {\cal L}_m} \epsilon \rho U^{\prime }}$ $k(TB\u2032\u2212Ts\u2032)d\u2032Lm\epsilon \rho U\u2032$ | 118.3 |

S | $\frac{\sigma \epsilon ^3}{\mu U^{\prime }}$ $\sigma \epsilon 3\mu U\u2032$ | 6.92 × 10^{−6} |

$\bar{K}$ $K\xaf$ | $\frac{kK}{{d^{\prime } {\cal L}_m}}$ $kKd\u2032Lm$ | 8.9 × 10^{3} |

h_{e} | $\frac{h^{\prime }_e}{d^{\prime }}$ $he\u2032d\u2032$ | 1 |

δ | $\frac{\alpha \mu U^{\prime }}{\epsilon ^2L(T^{\prime }_{\text{B}}-T^{\prime }_s)}$ $\alpha \mu U\u2032\epsilon 2L(TB\u2032\u2212Ts\u2032)$ | 4.66 |

A | $\frac{A^*}{{L^{\prime } d^{\prime }} \mu U^{\prime }}$ $A*L\u2032d\u2032\mu U\u2032$ | 2.14 × 10^{−6} |

$\tilde{k}$ $k\u0303$ | k/k_{c} | 1.17 |

P_{T} | $\frac{L^{\prime } \kappa _c}{U^{\prime } {{L^{\prime }_c}^2}}$ $L\u2032\kappa cU\u2032Lc\u20322$ | 0.76 |

Bi | $\frac{h_cd^{\prime }}{k}$ $hcd\u2032k$ | 0-0.1 |

### A. The tear film

#### 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 ⩽ *y* ⩽ *h*(*x*, *t*), we have conservation of mass, momentum (in the *x* and *y* directions) and energy, respectively:

The Reynolds number Re = ρ*U*′*L*′/μ ≈ 20, and together with ε = 10^{−3}, we have the reduced Reynolds number ε^{2}Re ≪ 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 ε^{2}RePr 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 −*L*^{′} ⩽ *x* < *X*^{′}(*t*^{′}), we assume that there is a thin fluid layer of constant thickness specified by

^{24,32,35}Here

*h*

_{e}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

At the film ends, we need boundary conditions as well. At those ends, the tear film is assumed to be pinned at a thickness *h*_{0}*d*; *h*_{0} 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 *h*_{e}*d* < *y*^{′} < *h*_{0}*d*, the tear film velocity components are equal to the motion of the ends. In the interval 0 < *y*^{′} < *h*_{e}*d*, 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:

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

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 forces^{52,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:

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

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.

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

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

*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, yielding^{16}

where

and

Here

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

### B. Heat transfer posterior to the tear film

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

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

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

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

Here *x*_{w} = 0.01 is chosen to be small, so that there is a narrow band for the transition between the two conditions.

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

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.

### C. End motion and boundary conditions for the film

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 Meuller^{76} 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 Δ*t*_{co} = 0.176; the downstroke, Δ*t*_{oc} = 0.082, and the interblink time, Δ*t*_{o} = 5. Thus, the blink cycle duration to be used in our experiment representing a normal eye blink is Δ*t*_{bc} = 5.258. A detailed expression of *X*(*t*) is given in Appendix A.

Braun and King-Smith^{31} explored full blink cycles using an ideal case: the sinusoidal lid motion

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

where *Q*_{top} and *Q*_{bot} 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

*Q*

_{top}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

Heryudono *et al.*^{32} explored the film dynamics’ dependence on *h*_{e}. We note that if *Q*_{top} = 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 *h*_{e} = 1 in our computations, corresponding to

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}

### D. The problem

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:

and

where *Q*_{bot} and *Q*_{top} depend on the choice of flux conditions. The initial condition for the film is specified as

Here *m* characterizes the shape of the polynomial and *h*_{min} is computed such that the volume *V*_{0} = *V*(0) under the curve is a specified value. The temperature field is initialized with body temperature:

### E. Numerical methods

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*) = −(*Sh*_{xx} + *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.

## III. RESULTS

### A. Sinusoidal lid motion with FPLM boundary conditions

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}, *V*_{0} = 1.576, λ = 0.1, *m* = 4, and *h*_{0} = 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 *h*_{e} = 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 *X*_{t}, but on the average the volume decreases over time due to evaporation.

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.

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-Smith^{31} 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

*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 *Q*_{c} 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 *Q*_{top} = −*X*_{t}*h*_{e}/2 + *Q*_{c} 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*(α + 2*k*π) − *V*(α)|, indicates that up to four digits of precision can be obtained with this method of evaporation compensation for several blink cycles.

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.

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 *Q*_{c} 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*, 2*t*_{bc} + *t*_{co}) − *h*(*x*, *t*_{bc} + *t*_{co})‖_{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}

### B. Realistic lid motion

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, *V*_{0} = 2.576, λ = 0.1, *m* = 4, and *h*_{0} = 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.

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.

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 *Q*_{c} 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*, 3*t*_{bc} + *t*_{co}) − *h*(*x*, 2*t*_{bc} + *t*_{co})‖_{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.

#### 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

*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 *Q*_{c} 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*, 2*t*_{bc} + *t*_{co}) − *h*(*x*, *t*_{bc} + *t*_{co})‖_{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.

We now turn to the dynamics of the OST, which is measured at the surface of the tear film.^{55,59} Li and Braun^{16} 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.

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}

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 drops^{59} 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.

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

^{73}

## IV. DISCUSSION

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 motion^{31} 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/cm^{2}/s in normals and 4.1 × 10^{−8} g/cm^{2}/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 experiments^{56,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).

## V. CONCLUSIONS

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}

## ACKNOWLEDGMENTS

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.

### APPENDIX A: END MOTION AND FLUX FUNCTIONS

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

Here Δ*t*_{bc} = Δ*t*_{co} + Δ*t*_{o} + Δ*t*_{oc} 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:

Here *Q*_{0lg} represents the supply of tear fluid from the lacrimal gland and *Q*_{0p} 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.)

Nondimensionalized parameters for lid motion and flux condition . | ||
---|---|---|

Parameter . | Description . | Value . |

Δt_{co} | Duration of upstroke | 0.176 |

Δt_{o} | Duration of interblink | 5 |

Δt_{oc} | Duration of downstroke | 0.88 |

Δt_{p} | Duration of punctal drainage | 0.5 |

Δt_{lg} | Duration of lacrimal gland supply | 0.088 |

t_{in} | Time location of influx peak | 0.176 |

t_{out} | Time location of outflux peak | 1.08 |

f_{top} | Fraction of supply/drainage from the superior side | 0.7 |

f_{bot} | Fraction of supply/drainage from the inferior side | 0.3 |

Q_{lg} | Tear supply from lacrimal gland | 0.337 |

Q_{p} | Drainage due to the puncta | 0.0297 |

Nondimensionalized parameters for lid motion and flux condition . | ||
---|---|---|

Parameter . | Description . | Value . |

Δt_{co} | Duration of upstroke | 0.176 |

Δt_{o} | Duration of interblink | 5 |

Δt_{oc} | Duration of downstroke | 0.88 |

Δt_{p} | Duration of punctal drainage | 0.5 |

Δt_{lg} | Duration of lacrimal gland supply | 0.088 |

t_{in} | Time location of influx peak | 0.176 |

t_{out} | Time location of outflux peak | 1.08 |

f_{top} | Fraction of supply/drainage from the superior side | 0.7 |

f_{bot} | Fraction of supply/drainage from the inferior side | 0.3 |

Q_{lg} | Tear supply from lacrimal gland | 0.337 |

Q_{p} | Drainage due to the puncta | 0.0297 |

### APPENDIX B: NUMERICAL METHODS

We define the pressure *p*(*x*, *t*) = −(*Sh*_{xx} + *Ah*^{−3}) as a new dependent variable. The thin film equations (14) and (15) then become

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

*x*and

*y*directions). At each time step, we used barycentric interpolation

^{87}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

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

The initial condition becomes

and the boundary conditions become

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.