{'fraction_non_alphanumeric': 0.03544437283965083, 'fraction_numerical': 0.020544066240943622, 'mean_word_length': 4.527037236913114, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computed Flow and Fluorescence Over the Ocular Surface

Fluorescein is perhaps the most commonly used substance to visualize tear film thickness and dynamics; better understanding of this process aids understanding of dry eye syndrome which afflicts millions of people. We study a mathematical model for tear film flow, evaporation, solutal transport and fluorescence over the exposed ocular surface during the interblink. Transport of the fluorescein ion by fluid flow in the tear film affects the intensity of fluorescence via changes in concentration and tear film thickness. Evaporation causes increased osmolarity and potential irritation over the ocular surface; it also alters fluorescein concentration and thus fluorescence. Using thinning rates from in vivo measurements together with thin film equations for flow and transport of multiple solutes, we compute dynamic results for tear film quantities of interest. We compare our computed intensity distributions with in vivo observations. A number of experimental features are recovered by the model.


Introduction
The purpose of the this paper is to construct and solve a model for the fluid flow and transport of solutes in the tear film on an eye-shaped domain. Using the computed solute distributions, we can calculate the intensity of light emitted from instilled fluorescein dye, and compare the computed result with a typical experimental observation. The results and comparison help understand tear film dynamics and imaging in vivo.
The tear film is critical for the eye's protection (Bron et al., 2004;Govindarajan and Gipson, 2010) and proper optical function (Montés-Micó et al., 2010;Tutt et al., 2000). The tear film is a thin liquid film with multiple layers that establishes itself rapidly after a blink. At the anterior interface with air is an oily lipid layer that primarily retards evaporation , which helps to retain a smooth well-functioning tear film (Norn, 1979). Posterior to the lipid layer is the aqueous layer, which consists mostly of water (Holly, 1973). At the ocular surface, there is a region with transmembrane mucins protruding from the epithelial cells of the cornea or conjunctiva. This forest of glycosolated mucins, called the glycocalyx, has been referred to as the mucus layer in past literature. It is generally agreed that the presence of the hydrophilic glycocalyx on the ocular surface prevents the tear film from dewetting (Tiffany, 1990a,b;Gipson, 2004). The overall thickness of the tear film is a few microns (King-Smith et al., 2004), while the average thickness of the lipid layer is on the order of tens to 100 nanometers (Norn, 1979;Yokoi et al., 1996;Goto and Tseng, 2003;King-Smith et al., 2011;Braun et al., 2015) and the thickness of the glycocalyx is a few tenths of a micron (Govindarajan and Gipson, 2010). This overall structure is reformed on the order of a second after each blink in a healthy tear film.
The aqueous part of tear fluid is primarily supplied from the lacrimal gland near the temporal canthus and the excess is drained through the puncta near the nasal canthus. Doane (1981) proposed the mechanism of tear drainage in vivo whereby tear fluid is drained into the canaliculi through the puncta during the opening interblink phase. Water lost from the tear film due to evaporation into air is an important process as well (Mishima and Maurice, 1961;Tomlinson et al., 2009;Kimball et al., 2010). We believe that this is the primary mechanism by which the osmolarity and other solute concentrations are increased in the tear film . Some water is supplied from the ocular epithelia via osmosis Cerretani and Radke, 2014;Braun et al., 2015).
The role of osmolarity on the ocular surface may be summarized as in Baudouin et al. (2013) and Li et al. (2016). The healthy tear film maintains homeostasis with the blood in the range 296-302 Osm/m 3 (or, equivalently, mOsm/L or mOsM) (Lemp et al., 2011;Tomlinson et al., 2006;Versura et al., 2010), while healthy blood is in the range 285-295 Osm/m 3 (Tietz, 1995). In dry eye syndrome (DES), the lacrimal system may be unable to maintain this homeostasis and osmolarity values in the meniscus rise to 316-360 Osm/m 3 (Tomlinson et al., 2006;Gilbard et al., 1978;Sullivan et al., 2010;Dartt and Willcox, 2013), and even higher values may be attained over the cornea. Using in vivo experiment and sensory feedback, Liu et al. (2009) estimated peak values of 800-900 Osm/m 3 . Similar or higher values were computed from mathematical models of tear film break up in Braun et al. (2015) or Peng et al. (2014), or for models of the whole ocular surface (Li et al., 2016). The estimates from these models are significant because the osmolarity in TBU, or in most locations over the ocular surface, is higher than in the inferior meniscus, where the only convenient method used in the clinic measures the osmolarity (Lemp et al., 2011).
The ongoing supply and drainage of tear fluid affects the distribution and flow of the tear film. A number of methods have been used to visualize and/or measure tear film thickness and flow, including interferometry (Doane, 1989;King-Smith et al., 2004, 2009, optical coherence tomography (Wang et al., 2003), fluorescence imaging (Begley et al., 2013;King-Smith et al., 2013b) and many others.
Fluorescein and other dyes have been used in a variety of ways, including: estimation of tear drainage rates or turnover times (Webber and Jones, 1986); assessment of the condition of the ocular surface via staining of epithelial cells (Bron et al., 2015, e.g.,); to visualize overall tear film dynamics (Benedetto et al., 1986;Begley et al., 2013;King-Smith et al., 2013b;Li et al., 2014); estimation of first breakup times of the tear film (Norn, 1969); and the temporal progression of tear film breakup areas (Liu et al., 2006). There are different ways that fluorescence may be used to visualize the tear film. In the dilute limit, the fluorescein concentration is below the critical concentration, and the intensity of the fluorescence from the tear film is approximately proportional to its thickness. In the concentrated (or self-quenching) limit, the intensity drops as the tear film thins in response to evaporation, and the thickness is roughly proportional to the square root of the intensity for a spatially uniform (flat) tear film Braun et al., 2014). In what follows, we will use FL to abbreviate fluorescence imaging to FL imaging. Fluorescein transport is difficult to measure, but we aim to better understand that and FL imaging via the mathematical model we develop here.
A variety of mathematical models have incorporated various important effects of tear film dynamics as recently reviewed by . The most common assumptions for these models are a Newtonian 3 of 38 tear fluid and a flat cornea (Berger and Corrsin, 1974;Braun et al., 2012). Tear film models are often formulated on a one-dimensional (1D) domain oriented vertically through the center of the cornea with stationary ends corresponding to the eyelid margins. We refer to models on this kind of domain as 1D models. Surface tension, viscosity, gravity and evaporation are often incorporated into 1D models; wetting forces have been included as well. Zubkov et al. (2012) formulated and studied a mathematical model that describes the spatial distribution of tear film osmolarity that incorporates both fluid and solute (osmolarity) dynamics, evaporation, blinking and vertical saccadic eyelid motion. They found that both osmolarity was increased modestly in the black line region of thinning near the meniscus (Miller et al., 2002) and that measurements of the solute concentrations within the inferior meniscus need not reflect those elsewhere in the tear film. This model gave smaller increases in osmolarity than spatially-uniform models Braun et al., 2015) because of the higher evaporation rates and lack of diffusion along the tear film in those latter models. Larger increases in osmolarity were reported in the eye-shaped domain model studied by Li et al. (2016) because they used larger evaporation rates as may be seen experimentally .
Compartment models of the tear film have also been developed. Gaffney et al. (2010) developed a model for the tear film with separate compartments for the central tear film, menisci and other locations with an impermeable cornea (no osmosis). The model allows water and osmolarity to pass between compartments so that mass is conserved, but it does not directly take into account fluid flow. From the model, they could estimate osmolarities obtained in the different compartments. Cerretani and Radke (2014) extended the model to include osmosis from the substrate.
Of interest in this paper are dynamics of the fluid flow and solute transport on an eye-shaped domain. To our knowledge, Maki et al. (2010a) were the first to extend models of fluid dynamics in the tear film to a geometry that approximates the exposed ocular surface. Besides specifying the tear film thickness at the boundary; they specified either the pressure (Maki et al., 2010a) or the flux (Maki et al., 2010b). Their simulations recovered features seen in 1D models such as formation of the black line, and captured some experimental observations of the tear film dynamics around the lid margins. Maki et al. (2010b) simplified the in vivo lacrimal supply and drainage mechanisms and imposed a time-independent flux boundary condition; under some conditions, they recovered hydraulic connectivity as seen experimentally. Li et al. (2014) improved that model by adding evaporation and a wettable ocular surface, as well as a time-dependent flux boundary condition that approximated the in-and out-flow of the aqueous layer of the tear film. Li et al. (2016) formulated a model that included osmolarity transport and osmosis from the tear/eye interface. This model provided, to our knowledge, the first global osmolarity distribution over the entire exposed eye. The permeability of the ocular surface was treated as either constant over the whole surface, or as a space-dependent function with lower permeability over the cornea and higher permeability over the conjunctiva. The location of the highest osmolarity varied with the permeability, and for the more realistic case with variability, the highest values were near the black lines over the cornea. This model did not include local tear break up dynamics, so they are an idealized case; however, the values found in that model are similar to recent calculations for tear breakup (TBU) models (Peng et al., 2014;Braun et al., 2015).
A number of mathematical models have been developed for TBU. Multi-layer models driven by dewetting van der Waals forces to rupture (corresponding to TBU in the eye literature) have appeared in linear (Sharma and Ruckenstein, 1986b), nonlinear Ruckenstein, 1985, 1986a) and non-Newtonian forms (Zhang et al., 2003(Zhang et al., , 2004. All of the theories could give reasonable TBU time (TBUT) ranges. Related papers, which may apply to eyes via analogy with lung surfactants, include Matar et al. (2002). For a comprehensive review of related bilayer work, see Craster and Matar (2009).
More recently, the DEWS report (Lemp et al., 2007) has argued that TBU is driven, at least in many cases, by evaporation, which causes increased osmolarity (concentration of ions) in the tear film, which may lead to irritation, inflammation and damage to the ocular epithelium. Currently, the authors are unaware of any method to directly measure the osmolarity in regions of TBU due to the very small volumes of tear fluid involved, the rapid dynamics there, and the extreme sensitivity of the cornea. Theoretical efforts have responded to this situation by creating models that incorporate evaporation and osmosis, as well as other effects, into the mathematical models to better understand the dynamics of the process at a small scale.
In  and Braun et al. (2015), a model for an evaporating, spatially uniform film was studied. The model was a single ordinary differential equation for the tear film thickness that included evaporation from the tear/air interface at a constant rate and osmotic flow from the tear/cornea interface that was assumed to be a semi-permeable boundary that allows water but not solutes to pass. They found that the model predicted equilibration of the tear film thickness and that the osmolarity could become quite large as the tear film thinned. Similar conclusions were found from the model of , which included van der Waals forces that stopped thinning at the purported height of the glycocalix; this allowed the model to be used at zero permeability at the tear/cornea interface.
These models were extended to include a specified space-dependent evaporation profile that varied in space by Peng et al. (2014) and Braun et al. (2015). In both models, the local thinning caused by locally increased evaporation led to increased osmolarity in the TBU region, which could also be several times larger than the isotonic value. Braun et al. (2015) also included fluorescein transport in their model, and could compute intensity distributions as well as fluorescein concentration. They gave preliminary results that showed that FL imaging techniques for TBU may need to be interpreted with care because fluorescein transport could alter the appearance of TBU from FL imaging relative to the actual thickness distribution. A detailed study that used an improved scaling has been been submitted . A model of TBU that incorporated an evaporation rate that depended on a insoluble surfactant concentration was studied by Siddique and Braun (2015); they concluded that more physics was needed in that model to correctly capture TBU dynamics.
In this paper, we investigate a model for tear film flow, evaporation, osmolarity and fluorescein transport, osmosis and fluorescence on an eye-shaped domain. The emphasis in this paper is on understanding the transport of solutes inside the tear film, and how that transport affects FL imaging of the tear film. Despite not specifically building in TBU into the model, the model will show similar effects from thinning and flow particularly in the black lines over the cornea. The paper is organized as follows. We briefly describe a sample experiment using FL imaging. We then develop the mathematical model for the fluid flow and solute transport on an eye-shaped domain. We then move on to results, followed by discussion and conclusion. The details of the mathematical formulation appear in the appendices.

A Sample Experimental Result
We begin by briefly describing the fluorescence imaging method used for visualizing tear film dynamics King-Smith et al., 2013a,b) for the images shown here; this description relies heavily on that used in Li et al. (2014). Video recordings were made from normal and dry eye subjects for 60 s after instillation of 1 µl of 5% (sodium) fluorescein. Subjects were instructed to blink about 1 s after the start of the recording and try to hold their eyes open for the remainder of the recording. The subjects eyes were illuminated with blue light and a blocking interference filter was used to reduce the response to reflected illumination light. The horizontal illumination width was 15 mm, thus including Computed Flow and Fluorescence Over the Ocular Surface 5 of 38 the cornea and part of the conjunctiva. The research protocol was approved by an Institutional Review Board in accordance with the Declaration of Helsinki. Informed consent was obtained from each subject at study enrollment.
Tear film thickness was studied by using the self-quenching of fluorescein, i.e., the reduction of fluorescent efficiency with increasing fluorescein concentration . When diluted by about the assumed standard 7 µl of tears, fluorescein concentration in the tears will be about 0.625%. As the tear film thins from evaporation, the fluorescent intensity is reduced inversely proportional to the square of concentration and is therefore proportional to the square of tear thickness. Thus, tear thickness is proportional to the square root of fluorescent intensity . The method estimates relative changes in thickness given an initial concentration, and does not directly visualize flow. However, the lubrication model that we use computes the film thickness directly, so it is appropriate to use this method to compare the theory with this type of experiment. Comparisons for open eye shapes have been made in our previous models (Maki et al., 2010b;Li et al., 2014).
In Figure 1, we show some images from the video recording of one subject (a different one than in our previous comparisons). The images show four different times post-blink. The transparent cornea appears darker than the surrounding reflective area corresponding to the conjunctiva; in these images, the contrast is enhanced to clearly show the features of interest. As time increases, the overall fluorescent intensity decreases away from the menisci. This is due to evaporative loss of water with the corresponding increase in fluorescein concentration. The self-quenching effect causes the intensity to decrease as the fluorescein concentration rises. The menisci remain bright throughout the experiment as they thin little and the fluorescein concentration changes little there as a result. Other sources for FL imaging results include, among others: Liu et al. (2006), Begley et al. (2013), Yokoi and Georgiev (2013), and Su et al. (2014). Localized effects such as TBU can occur as well; we don't deal with those dynamics in this paper, but FL imaging of TBU has begun to be studied theoretically . In this work, we aim to capture the overall tear film image dynamics in the absence of localized TBU.

Formulation
In this section, we present a mathematical model that incorporates fluid dynamics and transport of two solutes into a tear film model on a 2D eye shaped domain as shown in Figure 2. In Figure 2, (u , v , w ) are the velocity components in the coordinate directions (x , y , z ); z is directed out of the page and primed variables are dimensional. The gravitational acceleration g is specified in the negative y direction; we will neglect it in our computations in this paper.

Eye-shaped domain
The boundary curves of the eye-shaped domain are approximated from a digital photo by four polynomials. Two are parabolas in x and two are ninth-degree polynomials in y , and C 4 continuity is imposed where they join (indicated by dots) (Maki et al., 2010a,b;Li et al., 2014). s is the arc length of the boundary starting at the joint of the nasal canthus and upper lid, and is traversed in the counterclockwise direction as s increases. The unit vectors tangential and normal to the boundary curves are given by t b and n b , respectively. z = h (x , y ,t ) denotes the free surface of tear film and t is the time.

Governing Equations
We assume that the tear fluid is incompressible and Newtonian with constant density ρ, viscosity µ, specific heat c p , and thermal conductivity k. The tear fluid supplied from the lacrimal gland is more complex than water with salt ions, containing proteins and other organic molecules (Dartt et al., 2005;Dartt, 2009). Because these large molecules are present, whole tears have been found to be weakly shear thinning (Tiffany, 1991). Similar levels of shear thinning were incorporated by Jossic et al. (2009) and though some details change, the qualitative features of the flow and thickness distribution are, for our purposes here, similar to the Newtonian case. More recently, meibum (secreted from the glands in the eyelids) has been found to introduce elasticity into the tear film when temperature is below 32 • C (Leiske et al., 2011). Note that viscoelastic surfactants can affect film stability in a model system Bhamla et al., 2014Bhamla et al., , 2016. However, it is not yet clear how these measurements affect overall tear film properties at typical in vivo operating temperatures that are often near 35 • C (Efron et al., 1989). Because our current knowledge indicates that the shear thinning of the tear film is not strong (Tiffany, 1991) and because viscoelastic surface effects appear to be small at in vivo temperatures (Leiske et al., 2012), we make the simplifying assumption of constant properties.
We also assume the ocular surface is flat due to the fact that the characteristic thickness of the human tear film is much less than the radius of curvature of the ocular globe (Berger and Corrsin, 1974;Braun et al., 2012). The governing equations for the tear film thickness h (x , y ,t ), the osmolarity c (x , y ,t ) and the fluorescein concentration f (x , y ,t ) are derived from the incompressible Navier-Stokes and convection-diffusion equations, respectively. The two equations are coupled by the osmotic flux, that is, the fluid that flows from the ocular surface at z = 0 into the tear film driven by the osmolarity difference between the tear film and the cornea (Levin and Verkman, 2004). In the model, water is lost to the air, treated as a passive gas outside the tear film (z > h ), and water is supplied due to the subsequent increased osmolarity via this osmotic flux. It is assumed that fluorescein does not induce osmotic flux of water (see Appendix, or Braun et al., 2014).
After nondimensionalization and simplification using lubrication theory (e.g. Jensen and Grotberg (1993); Braun (2012)), we arrive at a system of PDEs for the dimensionless variables h(x, y,t), c(x, y,t) and f (x, y,t): We define the pressure as the evaporative mass flux J is given by and the fluid flux Q across any cross-section of the film is given by The conjoining pressure in modeling evaporation plays an important role in the tear film model. It is meant to mimic the effect of the glycocalix, whose transmembrane mucins are strongly wet by water and we assume that they arrest the thinning of the tear film. A secondary benefit is that the model allows solutions to be computed past the initial tear film breakup because the tear film thickness never reaches zero in this model. These aspects of evaporation competing with conjoining pressure are discussed by Winter et al. (2010) in the context of eyes, but the idea was developed by Potash and Wayner (1972) and Moosman and Homsy (1980). More recent versions of the approach may be found in Morris (2001) and Ajaev and Homsy (2001). The nondimensional parameters that arise are defined and given values in the following section and in Table 1. The dimensional parameters used in those expression are given in Table 2. A detailed derivation of the governing equations (3.1) -(3.3) can be found in Appendix A.
Once the fluorescein concentration f is computed, the FL intensity I can be computed from the thickness h and the fluorescein concentration f via the nondimensional version of (C.3) (Webber and Jones, 1986;Nichols et al., 2012) Here φ = κd f cr is the nondimensional Naperian extinction coefficient. For fixed thickness h, I decreases linearly for small concentrations; this is the dilute regime. Sufficiently far above a critical concentration f cr = 0.2%, the intensity asymptotes to a quadratically decreasing function with increasing f ; this is the self-quenching regime . For a tear film thinning by evaporation, the situation is slightly more complicated. For a flat tear film, the product h f is constant, and so in the dilute regime, I is constant during thinning, while in the self-quenching regime, the quadratic decrease for increasing f still holds . For a deforming tear film, mass conservation is more complicated, and it is a main purpose of this paper to compute how images of the tear film appear subject to our model for flow and solute transport.

Parameter Descriptions
The nondimemsional parameters we use are shown in Table 1; the dimensional parameters are in Table 2.
We note that when the ocular surface is permeable, we have a region corresponding to the cornea that is about 4.6 times less permeable than the surrounding conjunctiva, with a smooth monotonic transition between them. Details of the permeability function may be found in Li et al. (2016).

Initial Conditions
The initial conditions are the same as those of Li et al. (2016) subject to the addition of fluorescein concentration. For convenience we give all the initial conditions here. The initial thickness h(x, y, 0) is a numerically smoothed version of the function where x 0 = 0.06 and dist (X, ∂ Ω ) is the distance between a point with position vector X and a point on the boundary ∂ Ω (Maki et al., 2010a,b). It specifies a dimensional initial volume of about 1.805µl. This value is within the experimental measurements of 2.23 ± 2.5µl by Mathers and Daley (1996). The initial pressure p(x, y, 0) is calculated from equation (A.18) (Li et al., 2014). For the initial osmolarity, we assume the salt ions are well mixed and isotonic (302 Osm/m 3 , or 1 dimensionlessly), so that In addition, we specify the initial fluorescence concentration to be a constant, i.e.

Boundary Conditions
The boundary conditions follow those of Li et al. (2016) with the addition of no flux for fluorescein at the boundary. Along the boundary (denoted by ∂ Ω ), we prescribe constant tear film thickness: We set h 0 = 13, which is in the range of experimental measurement from Golding et al. (1997). In addition, we specify the fluid flux normal to the boundary via The details of this time dependent boundary condition are given in Appendix B and Li et al. (2016). For the osmolarity c(x, y,t), we employ a Dirichlet boundary condition We use a homogeneous Neumann boundary condition for f (x, y,t)  Here K corresponds to a nominal thinning rate of 4µm/min; to obtain a 20µm/min thinning rate, this quantity is reduced by a factor of 5.

Numerical Methods
To solve the problem numerically, we introduce the pressure p(x, y,t) as a dependent variable and substitute for J to arrive at (A.17)-(A.20). The corresponding boundary conditions must be applied, (3.11)-(3.14). Note that the flux condition (3.12) is readily converted into a Neumann condition for p. The initial conditions must be applied as well, using smoothed versions of (3.8) and (3.4), as well as (3.9) and (3.10). We use the Overture computational framework (http://www.overtureframework.org. Contact: W. D. Henshaw, henshw@rpi.edu), which is a collection of C++ libraries for solving PDEs on complex domains (Chesshire and Henshaw, 1990;Henshaw, 2002). The spatial discretization is carried out on five overlapping component grids, and the solution is found on the composite grid; further details are given elsewhere (Li et al., 2016).
To solve the equations (3.1)-(3.4), we first discretize the spatial derivatives using the second-order accurate finite difference method for curvilinear and Cartesian grids from Overture. Since the model equations (3.1)-(3.4) are weakly coupled via the solutes, the hybrid time stepping scheme implemented by Li et al. (2016) works well on this problem. To solve the coupled system, equations (3.2) and (3.3) are solved using an RKC method (Sommeijer et al., 1997). Then (3.1) & (3.4) are updated with c and f values and solved using a variable step size BDF method with fixed leading coefficient based on Brenan et al. (1996) and Maki et al. (2010a,b). The resulting nonlinear system of the BDF method is solved using Newton's iteration method. Further details are given in Li et al. (2016).

Results
We begin with an impermeable ocular surface and a nominal thinning rate of 4µm/min, which is at the upper end of normal rates (Nichols et al., 2005). We then proceed to a permeable ocular surface with either 4µm/min or 20µm/min nominal thinning rates, with different values for corneal and conjunctival permeability.
5.1 Impermeable ocular surface (P c = 0) In this subsection, we study the dynamics where the ocular surface is impermeable to water and solutes, P c = 0, and water evaporates via a uniform nominal (flat-film) thinning rate of 4µm/min. This case will result in the highest solute concentrations that we could see in this model because no water can be supplied by osmosis. Figure 3 shows the dependent variables h, c and f at times t = 5, 10, 15. We begin by discussing the change in the thickness h, corresponding to the first column. The initial thickness (not shown; see Li et al., 2014Li et al., , 2016) has a roughly flat interior with smoothly and rapidly increasing menisci at the edges of the domains (corresponding to the lid margins). The tear film menisci have a nondimensional thickness of 3 or more (15 µm or more dimensionally) in the pink band adjacent to the lid margins. Capillary-driven flow toward the menisci forms the so-called black line, as has been seen in many mathematical models (e.g. Miller et al., 2002;Braun and Fitt, 2003). In this case, the thinning from evaporation cooperates with capillarity to thin this region a little faster than it may otherwise occur. The interior of the tear film thins roughly at the specified thinning rate, while the thickness of the menisci may increase in places due to flow into that region inside the tear film. Figure 3 shows the change in the osmolarity c and fluorescein concentration f in the second and third columns, respectively. In both cases, the solute concentrations increase in the interior and black line regions, and change relatively little in the mensci. A large relative change may be expected from mass conservation if the relative change in the thickness is large from evaporation and the relative change is larger away from the menisci. For the osmolarity, the highest values are found near the black line regions, with a maximum value of c max = 1.83 for t = 15, similar to results found previously for osmolarity alone (Li et al., 2016). Thus, for the impermeable case, the maximum osmolarity is increased more than 80% from the isotonic value, to about 550 mOsM. For f , the distribution at each time is similar to c, but the maximum value obtained is larger, with f max = 2.09 being more than double the initial uniform concentration value. The maximum in f is larger because its diffusivity is about four times smaller than for the osmolarity. The larger diffusivity for osmolarity allowed more diffusion out of the high concentration regions into the surrounding tear film, thus lowering the peak value.
With the computed distributions of h and f , we can now use (3.7) to compute the intensity distribution from the fluorescein concentration. Recall that the FL intensity may be the most widely used imaging technique but that the intensity observed depends on both h and f . This implies dependence on both the initial thickness and f 0 , and the resulting intensity may then be more complicated than may be expected from commonly used interpretations Braun et al., 2014). Figure 4 shows the intensity for two different initial values of f . The smaller value is in the so-called dilute regime, and the larger is the self-quenching regime (Webber and Jones, 1986;Nichols et al., 2012).
The column on the left begins in the dilute regime where f < 1 (Webber and Jones, 1986;Nichols et al., 2012). The black line mentioned previously mentioned can be seen forming by t = 5, and this kind of imaging is how it got its name. The black line becomes that way because tangential flow into the menisci and the outer edge of the central region causing thinning and a subsequent reduction in intensity according to (3.7). A slight increase in I is seen on either side of the black line from the influx of fluid from the black line. As time increases, water is lost due to evaporation, f increases, but the intensity in the central region inside the black line changes relatively little. This is because the product f h that appears in the exponent is the most important term affecting the intensity in the dilute regime, and this product is nearly constant during evaporative thinning . A roughly unchanged I for the central tear film continues throughout this case because any changes in thickness and concentration are driven primarily by evaporation in this model.
Turning to the self-quenching case in the right column of Figure 4, we again see the black line develop by t = 5, with slight lightening on each side. However, the central region becomes a bit darker as time progresses. This is a feature of the self-quenching regime, where the intensity is dominated by the denominator of (3.7), and in evaporative thinning, f increases such that f ∝ 1/h and so f will substantiantially affect the intensity observed . For this case, the thickness does not change too much, but we will see more obvious differences in the subsequent cases.
There is an approximation of the intensity-thickness relationship that can be used fruitfully in the self-quenching regime. As suggested by Nichols et al. (2012), (3.7) can be approximated by where a is a constant typically chosen to calibrate the observed intensity to a known thickness for a given optical system Nichols et al. (2012). In our case, we choose a so that a √ I matches the initial thickness of unity in the central tear film region. Using our results for I and h, we tested how well this approximation worked for this simulation for two different initial values of f . The results are shown in Figure 5. If the approximation is working well, the difference |h − a √ I| should be close to zero. We see that for the dilute case on the left, each time has a large difference all over the domain. This is not unexpected for the dilute case, since the approximation is designed for the self-quenching regime. In that regime, the difference remains relatively small in the black line and central regions. The approximation did not work well in the menisci in either the dilute or self-quenching regimes, but this is to be expected because a was chosen to work best in the central and black line regions. We now proceed to the case where the ocular surface is permeable and the cornea is about four times less permeable than the surrounding conjunctival area Li et al. (2016). We consider two different evaporation rates: 4 µm/min and 20 µm/min. The former corresponds to the high end of the normal evaporation rates, while the latter is near the highest end of measured thinning rates King-Smith et al. (2010).

Evaporation rate 4 µm/min
The dynamics for h, c and f are shown in Fig.6. In the first column, h is seen to thin overall as time increases, and this is more noticeable away from the menisci. The black lines are still seen. The tear film becomes thinner over the cornea due to its lower permeability, and this appears as the circular central disk region as seen in our previous work without fluorescein (Li et al., 2016).
In the second and third columns, c and f increase as the tear film thins due to mass conservation of solutes. The osmolarity increases more over the cornea than the conjunctiva, and the osmolarity is highest in black lines over the cornea. The largest value appears over the cornea because of the lower permeability there (Li et al., 2016). The distribution of f is very similar in its pattern to c, despite the difference in diffusivity. In this simulation lasting 15 s, the maximum value is about 1.3 times the initial value, corresponding to about 390 mOsM peak concentration for the osmolarity or 1.3 f 0 for f . At this evaporation rate and with the supply of water via osmosis, the solute concentrations don't build up as much as the impermeable case, and because of this the diffusivities are insufficiently different to separate the solute distributions for this evaporation rate and end time. Fig. 7 shows the distribution of fluorescent intensity over the exposed ocular surface after using (3.7) together with the computed f and h. In either case, the meniscus stays at essentially constant intensity because h and f change little. The dilute case on the left has nearly no change in the central region inside the black line, but previously seen changes in the width of the meniscus is made visible by the fluorescence (Maki et al., 2010b;Li and Braun, 2012). The central region is similar to the impermeable case, and even though osmosis occurs, there is little tangential flow. As a result, the product h f remains nearly constant from mass conservation, so that I remains nearly constant Braun et al. (2014). For the self-quenching case (right), besides the black line, some overall dimming in the central region can be seen. Enhanced dimming is seen over the cornea, where the permeability is lower and the thinning is greater (see previous figure). Experimentally, we would expect the difference in dimming over the cornea compared to that over the conjunctiva to be unobservable due to the refective nature of the stroma beneath the conjunctiva. The drop in intensity over the central region is because f is larger than unity and the dominance of the denominator in Eq. 3.7 causes the intensity to drop with inverse square power of f Braun et al., 2014).
Despite the similarities of the solute plots, the intensity plots had some differences. These differences are highlighted by the difference between the computed thickness h and the approximation to the thickness from a √ I. Results for the two different initial fluorescein concentrations are shown in Figure 8. In either the dilute case with f 0 = 1/3 (left) or the quenching regime with f 0 = 2, the menisci have a large difference |h − a √ I| because a was chosen to work in the initially flat region outside the mensici. In the dilute case (left), the difference between computed and estimated values is large as well. In the quenching regime, the difference remains small in almost all of the black line and central regions. Near the nasal canthus where the black lines meet is the exception where the error becomes larger; this is not an area of measurement interest, to our knowledge. Thus the approximation to the thickness based on the intensity appears to be a good estimate of the thickness where the flow in the central region is small Nichols et al. (2012). 2 Evaporation rate 20 µm/min We now consider the variable permeability ocular surface with a nominal 20 µm/min thinning rate for the evaporation. This is a high evaporation rate, near the upper end of observed thinning rates in vivo . The dynamics for h, c and f are shown in Fig.9. In the first column, h is seen to thin rapidly outside of the menisci as time increases. The black lines are seen, but the central region steadily becomes thin as well. The tear film becomes even thinner over the cornea due to its lower permeability, appearing as the darker circular central region.
The thinnest tear film occurs in the inferior meniscus over the cornea.
In the second and third columns, both c and f increase outside of the menisci as the tear film thins due to mass conservation of solutes. The osmolarity increases significantly more over the cornea than the conjunctiva due to the lower corneal permeability to water, and the osmolarity is highest in the black lines over the cornea as in the lower evaporation rate case. However, the max value at t = 15 of 5.43, or 1630 mOsM, is much larger in this case. The distribution of f is similar in its pattern to c, with steadily increasing f outside the menisci and the highest values occurring over the cornea. However, the features in the f distribution are narrower than those in the corresponding c distribution because of the smaller diffusivity for f . The maximum f is higher as well, reaching 6.59 times the initial value f 0 . At this high evaporation rate, the simulation time is sufficiently long to show the distinction that develops between the solutes due to differing diffusion rates. Fig. 10 shows the distribution of fluorescent intensity over the exposed ocular surface after using (3.7) together with the computed f and h. In either case, the meniscus stays at essentially constant intensity because h and f change little. At this high evaporation rate, both the dilute (left) and selfquenching (right) cases show overall dimming, but it is much more dramatic in the latter case. In both cases, the dimming of the intensity is most over the cornea where the permeability is lowest and the tear film thins the most. The dilute case (left) starts out there, but f increases sufficiently, for the initial condition we used, to also enter the self-quenching regime later in the computation . The self-quenching case begins there, and the entire black line and central cornea become quite dark. It is not uncommon in quantitative experiments to enhance the FL images for processing so that better results can be obtained from images that are difficult to distinguish with the naked eye. The previously seen changes in the width of the meniscus is again made visible by the fluorescence (Maki et al., 2010b;Li and Braun, 2012).
Results for the two different initial fluorescein concentrations are shown in Figure 11 for |h − a √ I| with the high evaporation rate. In either the dilute case (left) or the quenching regime (right), the menisci have a large difference |h − a √ I| because a was chosen to work in the initially flat region outside the mensici. In the dilute case (left), the difference between computed and estimated values at essentially any location is large as well. In the self-quenching side, the difference between the two is smallest at t = 5 but remains relatively small particularly in the central region over the cornea. Thus the approximation to the thickness appears to be a good estimate of the thickness over the cornea in the high evaporation rate case Nichols et al. (2012).

Discussion
Our computations can capture a number of observed phenomena related to fluorescence imaging and tear film dynamics. The black line was named from its appearance in FL imaging studies, and the model captures its appearance. Many previous tear film models have captured the phenomenon from a fluid dynamics perspective, but we believe that we are the first to do so with a mathematical model combining fluid flow and fluorescence. In the menisci, the image remains bright because the tear film thins little and the FL concentration changes little. This seems to agree well with a number of experiments King-Smith et al., 2013a,b).
In the central region inside the black line, the tear film has relatively little tangential flow, and thins primarily by evaporation. The appearance of the thinning depends on the concentration of fluorescein and the rate of thinning. In the dilute regime, the evaporation causes the concentration to increase, and because the tangential flow is relatively small, mass conservation lead to h f being nearly constant. Then, the intensity changes little during thinning in dilute regime because the numerator of (3.7) is dominant . In the self-quenching regime, thinning again leads to increasing f , and while h f is nearly constant, the denominator of (3.7) is dominant so that the intensity decreases during thinning roughly proportional to f −2 . When the evaporation rate is high, the thinning may be dramatic, resulting in large changes in f , so much so that one may possibly switch from the dilute to the self-quenching regime during an experiment. Generally in that case, the self-quenching regime can be reliably used to estimate the thinning of the tear film from the intensity measurement.
From comparing the estimate for the thickness from a √ I with the computed thickness h, that the estimate worked best in the self-quenching regime ((where the approximation was derived by Nichols et al., 2012). Our work shows that this works well in the central region over the cornea for both slower and faster evaporaton. For spot measurement in the central region and fitting of the intensity measurements, it is possible for this match between thickness change and intensity to be very good . For quantitative measurement of thickness change in localized TBU, the situation is less clear if the initial fluorescein concentration is unknown (and it usually is not). Depending on the spot size of elevated evaporation, the relative rate of evaporation in the spot and initial fluorescein concentration, it is even possible for the intensity to increase for a time in the TBU region . The complication is that the tangential flow near TBU can rearrange the distribution of f because it doesn't diffuse as well as some other solutes in the tear film. However, dimming does occur when in the selfquenching regime, and so the approximation to determine the thickness from the intensity works at least qualitatively for TBU .

Conclusion
The mathematical model in this paper combines tear film flow, evaporation, transport of solutes (osmolarity and fluorescein), and osmosis on an eye-shaped domain representing the exposed ocular surface. We thus added fluorescein transport to our previous model (Li et al., 2016). Most of the increase in osmolarity still occurs over the cornea and in the black lines there. The results show that the location and value of the minimum tear film thickness and maximum solute concentrations are found to be sensitive to the permeability at the tear/eye interface.
There are some limitations to the current study for comparison with in vivo experiment; they are similar to those list for linking the osmolarity measured in vivo with the computed results in Li et al. (2016). We repeat those items here and add to them. (i) There is no lid motion to mix the tear fluid as occurs in vivo, which contributes to the variability observed in vivo. (ii) Once the black line is formed, there is little exchange between the tear film in the interior with the meniscus. The supply and drainage of tear fluid occurs in and affects primarily the meniscus, and together with a small relative thickness change, the meniscal values of the fluorescein concentration stay close to the initial value. (iii) The volume of tear fluid is probably large compared to aqueous-deficient DES subjects, and could affect the value of fluorescein concentration obtained, and thus the image observed. (iv) The initial value of the fluorescein concentration is not usually known for in vivo experiments. This quantity is clearly important for quantitative interpretation of fluorescence experiments (Webber and Jones, 1986;Nichols et al., 2012;Braun et al., 2014). (v) Similarly, for most experiments using fluorescein to measure relative thickness changes, the initial thickness is not known quantitatively. Knowing this quantity experimentally would greatly aid interpretation. (vi) The model does not include tear break up (TBU) per se, and there have been results that suggest that the osmolarity could be quite high in these localized regions (e.g., Liu et al., 2009;Peng et al., 2014;Braun et al., 2015Braun et al., , 2016. (vii) There is not staining of the ocular surface in this model, but it certainly can happen in the clinic (Efron, 2013;Bron et al., 2015) These points suggest directions for future research. Additional directions would include an dynamic lipid layer that affects the evaporation rate; possible models are those of of Bruna and Breward (2014) and Peng et al. (2014), as well as Stapf et al. (2016). Further directions would include blinking, where mixing of the tear fluid could occur in the model; this would allow closer comparison for both osmolarity measurement and imaging. Localized TBU in the model would also be desirable, but being able to resolve small TBU areas in the domain would almost certainly require adaptive mesh refinement (Henshaw and Schwendeman, 2008). the permeability is the same for these ions as naturally occuring components of osmolarity, but they possess different diffusivities ( Table 2).
Pe c = U 0 L D c is the Pèclet number for salt, P c = P tiss v w c 0 εU 0 is the nondimensional permeability of the ocular surface and B = 3 f cr /c 0 is the ratio of osmotic components. The tissue permeability P tiss will take on different values as described in Section 5. We now examine the relative sizes of the osmotic suppy terms in W o . The critical concentration of fluorescein is 0.2% in water; this is a mass fraction of 0.002 g FL/g water. We assume that the aqueous component of the tear film has the same properties as water. There are roughly 1000 g/l of water, so that f cr = 2 g /liter of water. The molar mass of Na 2 FL is 376 g/mole, so the concentration of the undissociated molecule f cr = 0.0053 mole/liter or 5.3 millimolar (mM). The concentration of ions is then three times the concentration of FL 2− , hence the factor of three in B. Substitution yields B = 0.0525. Because B is so small and f is at most only a few times the size of c, we have neglected the contribution of f to W o . Thus, we will use W o = P c (c − 1) in this paper.
Applying lubrication theory by neglecting all the small terms for the fluid equations, we then have the following leading order approximations for the water, momentum and energy conservation.
In 0 < z < h(x, y,t): For the osmolarity, we expandc(x, y, z,t) as The leading order equation for the osmolarity is which impliesc 0 is independent of z, i.e.c 0 =c 0 (x, y,t). We proceed to the next order so as to find an equation forc 0 , and we obtain The equation resulting fromf is similar. We use ∇ = (∂ x , ∂ y ) and ∆ = (∂ 2 x + ∂ 2 y ) to represent the differential operators applied on them for convenience since h andc 0 are all independent on z. We note that since p is independent of z, the leading order value throughout the film is specified by (A.11).
At z = 0, using the simplified expression for W o yields: For the tear film, we solve for the velocity and temperature fields, integrate the mass conservation equation and use the kinematic condition to derive a PDE for h(x, y,t): For the osmolarity, we integrate Equation (A.7) with respect to z from 0 to h. Noting thatc 0 = c 0 (x, y,t), we then have ∂ zc1 (x, y, h,t) − ∂ zc1 (x, y, 0,t) = Pe c [h∂ tc0 + ∇c 0 · Q] − h∆c 0 . Parameters in the formulations are tabulated in Table 3. On-ramp location for lacrimal gland peak 4.2 s lg,o f f Off-ramp location for lacrimal gland peak 4.6 ∆ s lg On-ramp and off-ramp width of lacrimal peak 0.2 p out Fraction of drainage from upper punctum 0.5 s p,lo Lower punctal drainage peak location 11.16 s p,up Upper punctal drainage peak location 11.76 ∆ s p Punctal drainage peak width 0.05

C. Fluorescence Model
Fluorescein, when discussing the tear film, is typically referring to the use of the salt Na 2 FL, which ionizes to 2Na + and FL 2− at physiological pH values of 7.2-7.6 for tears (Doughty, 2010;Efron, 2013;Bron et al., 2015). It glows green when exposed to blue light (Webber and Jones, 1986), and has been used in a variety of ways for the tear (Efron, 2013;Bron et al., 2015), including as a first indication of TBU (Norn, 1969), as well as a variety of biological applications (Lakowicz, 2006). The fluorescent intensity observed in the tear film is a product of the efficiency of fluorescence E f ( f ) and the absorptance A f (h , f ) . Here h is the tear film thickness and f is the fluorescein concentration. For the efficiency, E( f ) = C 1 1 + ( f / f cr ) 2 , (C.1)