## Abstract

A new modeling approach is proposed for representing a subpolar ocean, whose upper layer is partly convected with the lower layer. A simple box model, built as a reference, has one active box to be the upper layer of a colder ocean, which interacts with the warmer box and the lower box. The active box receives atmospheric forcing (cooling and precipitation) and a parameterized freshwater (or sea ice) flux as well, while the other boxes have their properties fixed. The active box, interacting with the warmer box, possesses a thermal-driven state, at which the warmer water enters the active box, is cooled by the atmosphere, and becomes denser. The lower box adds another solution: a convected state appears in the vicinity of the nonconvected state. The nonconvected state either separates from or is absorbed into the convected state; that is, the entire upper box is convected with the lower box, once the lower-box density is close to the upper-box density.

The new component to the simple model is a probability distribution function (PDF) on the temperature–salinity (*T*–*S*) plane for the active box. The PDF in this probability box model represents heterogeneity in the upper layer, whereas one box has to be homogeneous in an ordinary box model. A *T*–*S* distribution retains only the probabilities of different water types, while their locations are discarded. The mechanisms to increase and reduce heterogeneity are modeled by the divergence and convergence of the PDF on the *T*–*S* plane, respectively. The heterogeneity is generated by the intrusion of exterior water as well as variability in the atmospheric forcing and freshwater flux, while the heterogeneity is reduced by horizontal diffusion within the box. Convection with the lower box tends to concentrate the PDF to the *T, S* of the lower box. Under the exterior condition that could produce both nonconvected and convected states in the simple box model, there is only one state of the upper box, which is partly convected, in the probability model. This intermediate state is possible when the divergent mechanism is intense, and the convergent mechanism is weak. Thus, the on–off convection in the simple box model is replaced with an intermediate state between the convected and the nonconvected states. It is suggested that, once mesoscale variability maintains heterogeneity, convection in the subpolar ocean is more robust against freshwater input.

## 1. Introduction

It has been suggested that the present global deep circulation has the North Atlantic Deep Water (NADW) as a major source in the northern North Atlantic (Broecker 1990). The NADW is formed from the warm and saline Atlantic water (AW), which flows northward and descends after cooling as a typical thermal-driven thermohaline circulation. The NADW is significantly contributed by the Denmark Strait Overflow Water (DSOW) (Luyten et al. 1993), which flows out of the Greenland, Iceland, and Norwegian (GIN) Seas. In the GIN, as schematically shown in Fig. 1, the DSOW is the product of the Greenland Sea Deep Water (GSDW) and the Norwegian Sea Deep Water (NSDW). The GSDW is the heaviest water in this region, being directly formed by deep convection to the ocean bottom (Aagaard and Carmack 1989). In consequence of thermohaline circulation, the downward water transport is partly responsible for the properties of the GSDW and the NSDW.

However, either this thermohaline circulation in the northern North Atlantic or the deep convection in the GIN has not been persistent in the past, as indicated by paleoclimatic records. The alternative state of the present thermohaline circulation is a salinity-driven state, at which fresh- and cold water flows southward near the sea surface. The most recent such event occurred during the Younger Dryas, when freshwater from deglaciation poured into the North Atlantic (Lehman and Keigwin 1992). Another freshening event was observed to a much less extent in the GIN and the northern North Atlantic; that is, during the period of the Great Salinity Anomaly (GSA), the GIN was covered by low salinity water, and deep convection ceased after 1968 (Dickson et al. 1988). This was the period when a nonconvected state occurred as the alternative of the convected state. Chemical tracer data suggest that deep-water formation in the GIN has been reduced greatly since the GSA (Rhein 1991; Schlosser et al. 1991). The DSOW may be reduced after 1968, and the NADW would be modified gradually. Hence, a clarification of the responses of the GIN to a freshwater input will be a forward step to the estimation of variability in the global deep circulation.

Similar switching between the convected and the nonconvected states was observed in the Labrador Sea. Hydrographic data collected at OWS Bravo (Lazier 1980) showed that no convection reached below the 100-m depth during the GSA (1968–71), but the water column was otherwise convected to 400-m to 1500-m depths. Thus, the deep convection ceased for several years once the Labrador Sea was covered by the GSA, and then it recovered afterward.

A simple two-box model proposed by Stommel (1961) has shown that the thermohaline circulation can possess a bimodal solution consisting of the thermal- and the salinity-driven states. There are two mechanisms essential to the bimodal solution: one is a difference in timescales between (quick) thermal response and (slow) salinity response to the atmosphere and the other is a nonlinear exchange rate between the boxes as a function of the density difference. The switching between the two states occurs after a sudden input of freshwater or salt into one of the boxes.

A bimodal solution exists also between the convected and the nonconvected states. It is simply shown that the bimodality is possible, when the upper ocean cools and receives freshwater such as a subpolar ocean (Welander 1982; Ikeda 1987). Once convection ceases, the upper ocean continues to be freshened without mixture of a saline water from the lower ocean. In contrast, when convection mixes the upper ocean with the lower ocean, cooling increases the upper-ocean density and maintains convection continuously.

Here we have the following open question: whether the GIN will return to the convected state. Welander (1982) proposed an oscillatory behavior between the convected and the nonconvected states using a vertical two-box model. In addition to the timescale difference between the thermal and the salinity responses to the atmosphere, another key condition is that the atmosphere tends to increase the temperature and salinity of the upper layer. This condition is possible only in the tropical regions, whereas the oscillation does not occur in a subpolar ocean. In the present work, a focus is given only to the convected–nonconvected bimodality because it seems more sensitive to a sudden freshwater or salt input than the thermohaline bimodality.

Adding a stochastic component to Stommel’s (1961) box model, Cessi (1994) has suggested that a fluctuating external force makes one of the two possible states more probable than the other. Cessi (1996) has shown that Welander’s (1982) box model also has an oscillatory solution responding to a fluctuating force even in the nonoscillatory range with a steady force. Thus, her models indicate that a simple box model, having only deterministic solutions, might provide an unrealistic solution.

Spatial variability is expected to introduce effects similar to temporal fluctuations. Although it is straightforward to include time-dependent, fluctuating components in a box model, the spatial variability might seem contrary to a box model. In the present paper, this heterogeneity is represented by a probability distribution function (PDF) defined on a temperature–salinity (*T*–*S*) plane; that is, a single set of *T, S* in the ordinary box model is extended to the PDF on the *T*–*S* plane. From the viewpoint of three-dimensional modeling, a *T, S* value at each location is dropped, whereas only a rate of the water with a certain *T, S* value is retained. The mechanisms that increase or reduce the heterogeneity are explicitly included in a 3D model, whereas they are parameterized by the divergence or convergence of the PDF in the present model.

An ordinary box model is introduced in section 2 for the upper layer of a colder ocean by one active box that cools and interacts with the warmer box and the lower box. The box model results are shown, producing the convected and the nonconvected states. On the basis of the simple box model, the probability model is developed in section 3, and its results are presented in section 4. The model behaviors are discussed in section 5.

## 2. Simple box model

### a. Model

A box model is first generated as a basis for a probability model. The configuration of this box model is shown in Fig. 2. Here the model represents a subpolar ocean in general, while the GIN is kept in mind as an example. The model has only one active box of the upper layer in the subpolar ocean, which interacts with the two domains that represent the ambient lower ocean under the active box as well as the warmer domain adjacent to the active box. The model with only one active box may be justified by the objective in the present paper: a role of spatial variability is examined on the bimodality of a convective ocean using the newly developed probability box model. The box also receives atmospheric cooling, precipitation (or evaporation), and freshwater flux from the polar region (or in the form of sea ice).

Temperature and salinity, *T* and *S,* in the box are prognostically calculated

where *t* is time. The subscripts *a, o,* and *w* denote the atmosphere, the lower ocean, and the warmer domain respectively. Hereafter, the freshwater (and sea ice) flux is included in the atmospheric forcing. Coefficients *k*_{T}, *k*_{S}, *k*_{o}, and *q* represent interaction rates with the respective external domains. The density, which is necessary to determine *k*_{o} and *q,* is calculated as a linear function of *T* and *S*

where only a density anomaly is considered with the coefficients of

Following Munk and Anderson (1948), the vertical interaction coefficient *k*_{o} is given the form

where *E* is chosen to be 2 × 10^{−10} kg^{1.5} m^{−4.5} s^{−1} and *ρ*_{m} is 0.001 kg m^{−3}. The same form was used by Ikeda (1987) for the Labrador Sea box model. The value chosen for *E* corresponds to a vertical diffusion coefficient *k*_{o}*H*^{2} = 2 × 10^{−4} m^{2} s^{−1} with a typical value of *ρ*_{o} − *ρ* = 0.4 kg m^{−3} and the upper-layer thickness of *H* = 500 m. The maximum value of *k*_{o} is set (∼3 × 10^{−5} s^{−1}) by (2.3), as *ρ* exceeds *ρ*_{o} − *ρ*_{m}; that is, the box is mixed with the lower ocean in a few days.

The volumetric exchange coefficient *q* has the form

which is proportional to the density difference between the box and the warmer domain, as suggested by Stommel (1961). This form is identical to Stommel’s, whereas there is a difference in the processes represented by them. Stommel (1961) considered the thermohaline circulation, in which warm water flows poleward in the upper ocean, descends after cooling, and flows equatorward in the lower ocean. Each of his two boxes included the entire water column. In the present study, a water exchange within the upper ocean is modeled. It is likely that a frontal current exists between the cold and the warm domains, inducing water exchange across the front due to mesoscale eddy formation (Griffiths 1995). Since the frontal current is more intense with the larger density difference, the exchange rate is reasonably assumed to increase as the density difference becomes larger. This parameterization has been used by Gargett and Ferron (1996).

The other model parameters are chosen for a typical subpolar ocean similar to the GIN (Fig. 2). However, it is noted that an important condition is the thermal-driven state appearing near the lower ocean property. Hence, only the orders of the magnitudes are concerned with here. The box includes an area of 2 × 10^{12} m^{2} (∼1200 km × 1500 km). The upper-layer thickness *H* is chosen to be 500 m (Aagaard et al. 1985; Clarke et al. 1990). The reference values of *T* and *S* are taken to be −1°C and 34.7 ppt, and the lower ocean property is set around the reference values. Hereafter, the temperature and salinity values are measured relative to the reference values.

The values of *T*_{a}, *S*_{a} and the interaction coefficients, *k*_{T} and *k*_{S}, are chosen based on estimates of air–sea interactions and the export of sea ice and freshwater from the polar ocean. The annual mean heat flux is estimated to be ∼−100 W m^{−2} (McCartney and Talley 1984), where a negative sign denotes flux from the ocean to the atmosphere. The flux is distributed over the 500-m layer, corresponding to *k _{T}*(

*T*) ∼ −5 × 10

_{a}− T^{−8}K s

^{−1}. Based on a realistic temperature difference of 5°C (Table 1),

*k*

_{T}is chosen to be 10

^{−8}s

^{−1}. Thus, the timescale of an oceanic response may be a few years: 1/

*k*

_{T}∼ 10

^{8}s.

As suggested in various papers (e.g., Bryan 1986), salinity feedback is much weaker than temperature feedback: *k*_{S} ≪ *k*_{T}. Hence, *k*_{S} is chosen to be 0.03*k*_{T} = 0.03 × 10^{−8} s^{−1}; that is, a salinity response occurs in 100 years. The salinity difference is taken to be −10 ppt (Table 1), yielding the salinity flux to be *k*_{S}(*S _{a}* −

*S*) ∼ −3 × 10

^{−9}ppt s

^{−1}. For the subpolar ocean, the freshwater flux due to precipitation minus evaporation is negligible (Baumgartner and Reichel 1975), whereas the export of sea ice and/or fresher water is a common feature. Here, the salinity flux is equivalent to 10

^{5}m

^{3}s

^{−1}of the freshwater export from the polar ocean, being similar to the estimate made by Aagaard and Carmack (1989) for the GIN Seas.

The water property in and the exchange rate with the warmer box are chosen as (*T*_{w}, *S*_{w}) = (5, 0.5) relative to (−1°C, 34.7 psu), and *C* = 3 × 10^{−8} kg^{−1} m^{3} s^{−1} (Table 1). Under the condition that the active box is at the thermal-driven state (section 2b), *q* becomes 8 × 10^{−9} s^{−1}. This value yields the volume exchange rate of 7 × 10^{6} m^{3} s^{−1} with the warmer box. The water property and the exchange rate are comparable with the warm and salty subpolar mode water through the Faroe–Shetland Channel and Denmark Strait into the Norwegian Sea (McCartney and Talley 1984; Aagaard and Carmack 1989).

### b. Results

In this section, model results are presented first without vertical mixing to confirm that the active box property becomes similar to the reference values. When the vertical mixing is excluded (*E* = 0), the solution should essentially be equivalent to the bimodal solution obtained by Stommel (1961). Starting with an arbitrary initial condition, a solution converges to either (*T, S*) = (−3.064, −0.666) or (−0.545, 0.123), where (*T, S*) are relative to (−1°C, 34.7 psu). These two states are referred to as equilibrium states 1 and 2 (ES-1 and ES-2): ES-1 is colder and fresher, and ES-2 is warmer and saltier. At ES-1, the box is lighter than the warmer box; that is, it is the salinity-driven state. The box is heavier at ES-2, or the thermal-driven state.

The evolution of the state in the box is displayed on the *T*–*S* plane in Fig. 3, starting with various initial conditions. In this paper, *T* and *S* are taken as the vertical and the horizontal axes, respectively, following the conventional way for physical oceanography. The convergence toward ES-1 or ES-2 slows down as the state approaches ES’s. The water exchange with the warmer area is much weaker at ES-1 than at ES-2 because |*ρ*_{w} − *ρ*| is much smaller at ES-1 (∼10% of that at ES-2). Hence, ES-1 is colder and fresher than ES-2. The bimodal character is fundamentally consistent with that found by Stommel (1961). In this paper, a focus is given only to the thermal-driven state, ES-2.

Sensitivity of the bimodal character is examined with *k*_{S} and *C* varied. It is shown that ES-1 exists always, while ES-2 disappears as the freshwater flux is increased, or as the exchange with the warmer domain is reduced: the common tendency is a reduction in *S.* When only one parameter is varied with the others kept constant, the transition from the bimodal to the unimodal solutions occur around *k*_{S} = 0.05*k*_{T} and *C* = 10^{−8}. This transition is qualitatively similar to that in Stommel (1961).

Vertical mixing between the upper and the lower boxes is now included: *E* is changed to 2 × 10^{−10} kg^{1.5} m^{−4.5} s^{−1}. The evolutions are shown in Fig. 4 for two cases of the lower-box property; that is, (*T*_{o}, *S*_{o}) = (0, 0.3) and (0, 0.5). The vertical mixing changes the evolution particularly in the subplane for *ρ* > *ρ*_{o}, where the solution rapidly converges to (*T*_{o}, *S*_{o}). In the subplane for *ρ* < *ρ*_{o} also, once the density of the box approaches the lower-layer density *ρ*_{o}, the solution quickly converges to (*T*_{o}, *S*_{o}). This state is to be distinguished from ES-2 and is referred to as ES-3 or the thermal-driven convected state.

The distinct difference between the two cases in Fig. 4 is whether ES-2 exists or not. Once the lower-box density is close to that in the upper box, convection occurs, and ES-2 is combined with ES-3. The critical value of *S*_{o} exists between 0.35 and 0.36, as long as the other parameters are fixed. Thus, the on–off convection is realized in the present box model, and will be examined in section 4 using the probability box model. It is noted that the convection shut-off may be related to an actual oceanographic event such as the GSA: the external forcing is modified to make the smaller upper- box density available, and a sudden reduction of the density is induced in the upper box, shutting off convection.

We now discuss whether the water exchange in the present model is appropriate for the subpolar ocean. Only vertical convection is included in the exchange between the active box and the lower box. The exchange between the active box and the warmer domain is represented by eddy-induced mixing or inflow/outflow in the upper ocean. The component missing in the model is the vertical overturning associated with thermohaline circulation, which consists of a poleward flow in the upper ocean, a descent in the subpolar ocean, and an equatorward flow in the lower ocean, at the thermal- driven state. The contribution of the poleward flow can be included in the model by parameterization of the volumetric exchange coefficient. The descent would carry the water in the active box into the lower box. Since the lower box property is fixed in the present model, no fundamental difference is produced by the thermohaline circulation.

## 3. Probability box model

The upper layer of the subpolar ocean, which usually has heterogeneity, has been treated as a box, filled with a homogeneous water mass in the simple box model (section 2). This spatial uniformity could be supported by uniform atmospheric forcing and uniform distribution of the exchanged water. However, in the real ocean, spatial variability of the atmospheric forcing exists, reflecting physical processes such as cooling due to a synoptic atmospheric low pressure system. Sea ice is exported from the adjacent shelf by synoptic winds and oceanic mesoscale eddies. The exchange with the warmer domain occurs on the boundary, injecting blobs of the warmer water, which wander in the subpolar ocean as sources of heterogeneity.

Horizontal variability is a common character in a convected area. Johannessen et al. (1991) have observed a chimney with a diameter of ∼10 km in the Greenland Sea. The chimney sat in a cyclonic eddy, implying that the eddy was a precondition for convection. The horizontal variability associated with mesoscale eddies are well exhibited in a numerical model with a high resolution (Heburn and Johnson 1995). Their horizontal scales are 10–100 km and timescales are weeks to months.

It is clear that a three-dimensional model is required to describe the structures and evolution of heterogeneity explicitly, although it is also obvious that a high-resolution three-dimensional model requires large computer resources. Here, we take a different approach: a probability model, in which the probability of each water property is introduced as a new variable to count heterogeneity. However, we have to drop information on the position of each water property. We also have to prescribe mechanisms of generation and dissipation of the heterogeneity. The probability model is built on the simple box model. Instead of single values of *T* and *S* for the box, a probability distribution function (PDF) on the *T*–*S* plane is introduced. Evolution of the PDF is determined from the external forcing and internal mixing.

As shown in Fig. 5, the PDF, *P*(*T, S*), is larger for higher probability. The mean values, [*T*] and [*S*], are integral of *T* and *S* weighted by *P,* respectively. A discrete form is used here by discretizing the *T*–*S* plane: *P* = *P*_{ij} at (*T*_{i}, *S*_{j}) ≡ (*i*Δ*T, j*Δ*S*), where (Δ*T,* Δ*S*) denotes a grid size on the *T*–*S* plane. The spatial mean is expressed as

where the sum of the probability becomes unity based on its definition:

The evolution of *P* is categorized into three kinds of mechanisms. The external forcing without variability mainly induces a shift of peak on the *T*–*S* plane. Variability in the external forcing and intrusion of the external water tend to spread *P* and reduce the peak amplitude, while internal mixing tends to narrow *P* and increase the peak amplitude. The formulation is first made for the evolution of *P*_{ij} due to the external forcing that generates no variability; that is, it is the evolution described in the simple box model. Before showing the form of *P*_{ij} evolution, we derive the time derivatives of *T* and *S* that would occur if the model were treated as the simple box model filled with homogeneous water

where *k*_{o} is determined from (2.3) with the local value of *ρ* as a function of (*T*_{i}, *S*_{j}). However, *q* is computed from (2.4) with *ρ* replaced with the mean [*ρ*]. The choice of [*ρ*] reflects the assumption that the horizontal exchange is dependent on a difference between the mean density in the box and that in the warmer domain. An alternative may be the individual density *ρ*_{ij}, although this choice is not taken here, because the horizontal exchange rate is expected to be controlled mainly by the large-scale density structures.

Those time derivatives should be accounted for by the time derivatives of *P*_{ij}. The evolution is schematically described in Fig. 6. An increase in *T,* for example, is equivalent to conversion of water from (*T*_{i}, *S*_{j}) to (*T*_{i+1}, *S _{j}*); that is, some part of

*P*

_{ij}is converted to P

_{i+1,j}. For both signs of the time derivatives of

*T*and

*S,*only the parts related to the conversion from

*P*

_{ij}are first formulated in the time derivatives of

*P*

_{ij}and its four neighboring grids,

It is noted that contributions of all grids to their four neighboring grids are summed up, and hence, the time derivative of *P*_{ij} has a contribution of *P*_{i−1,j} or *P*_{i+1,j}, as well as that of *P*_{i,j−1} or *P*_{i,j+1} in addition to itself. It is shown in the appendix that the evolution governed by (2.1) is recovered from (3.3).

The convection is treated separately from the formulation introduced above. Only for the grids with the density greater than a certain value close to the lower- box density, the convection mixes the upper layer instantaneously with the lower layer (Fig. 6):

where *ρ*_{m} = 0.001 kg m^{−3} and *P*_{IJ} denotes the PDF corresponding to the lower-box property. The choice of this critical density *ρ*_{o} − 4*ρ*_{m} is related with the mixing timescale determined from (2.3) used for the simple box model: *k*_{o} ∼ 10^{−6} s^{−1} at *ρ*_{o} − *ρ* = 4*ρ*_{m}, where *E* is chosen to be 2 × 10^{−10} kg^{1.5} m^{−4.5} s^{−1}. Thus, the mixing occurring with a timescale shorter than 10 days in the simple box model is replaced with the instantaneous mixing.

Among an evolution of heterogeneity, its generation is formulated first. The atmospheric forcing has variabilities with horizontal scales larger than a few 100 km and may not be effective in generating heterogeneity. However, sea ice and freshwater are expelled from the polar region partly due to oceanic mesoscale features. Hence, their contribution is more effective. The generation due to these mechanisms is represented by the PDF spreading over the *T*–*S* plane; that is, the transfer from a higher value to a lower value proportional to its gradient (Fig. 6). The components related to the gradient between *P*_{ij} and *P*_{i+1,j} and the gradient between *P*_{ij} and *P*_{i,j+1} are

It is noted that (3.4) is summed up for all grids, and the time derivative of *P*_{ij} is contributed by *P*_{ij} itself and the four neighboring girds. Once the coefficients *D*_{T} and *D*_{S} are set to be constant, a continuous form equivalent to (3.4) is represented by the second derivative terms, ∂^{2}*P*/∂*T*^{2} and ∂^{2}*P*/∂*S*^{2} (appendix).

The other generation mechanism, which is a horizontal exchange with the warmer domain, could be represented by adding P corresponding to its water property. However, this point is far from the mean property in the box and requires a large *T*–*S* domain (and many grids). Hence, only its essential role is taken here. A nonuniform property in the warmer domain generates heterogeneighty in the active box. In addition, nonuniform mixing (even with a homogeneous warmer domain) is a source of heterogeneighty, as the warmer water intrudes in the box interior but is not completely mixed. Similar to effects of atmospheric forcing and freshwater flux, the PDF spreads. These effects are combined into (3.4).

It is hard to evaluate the effects of the actual generation mechanisms. Hence, the approach taken in the present study is to estimate *D*_{T} and *D*_{S} to the best of our knowledge, and to examine the model sensitivity of *D*_{T} and *D*_{S}. Variabilities may be fractions of the external forcing terms, and responses of the subpolar ocean to the variabilities are further reduced, because oceanic responses work on the variabilities as a low-pass filter. The mean heat and salt forcing terms are shown in section 2 as

The magnitudes of the generation mechanisms for *T* and *S* are chosen to be 0.1% and 0.2% of the mean heat and salt forcing terms, respectively. The generation magnitudes are divided by the scales of ∂*P*/∂*T* and ∂*P*/∂*S,* yielding *D*_{T} and *D*_{S}. Once the ranges of the PDF are set 1°C and 0.5 ppt, the gradients become 1 K^{−1} and 2 ppt^{−1} with respect to temperature and salinity. Hence, the reference values are

The internal mixing in the upper layer is formulated as the last component. In the mixing process, two water types produce a new water type with a property between the original properties: *P*_{i′j′} and *P*_{i"j"} are transferred to *P*_{ij} at (*i, j*), which can be any grid between (*i*′, *j*′) and (*i*″, *j*″). However, in the present study, (*i, j*) is specified to be the midpoint for simplicity (Fig. 6), reflecting the assumption that the new water type is produced from equal amounts of the two original water types. In the following, it is specified which grids are mixed with each other.

The typical case considered with the probability model is that convected patches exist among a nearly homogeneous, nonconvected portion in the box. Hence, the water having (*T*_{o}, *S*_{o}) is assumed to contact with all the other types of water. The mixing between (*I, J*), which corresponds to (*T*_{o}, *S*_{o}), and the other grids (*i*′, *j*′) are

When *I* + *i*′ is an odd number (2*i* + 1), this term is divided by 2 and distributed to *dP*_{ij}/*dt* and *dP*_{i+1,j}/*dt* equally. If *J* + *j*′ is also an odd number (2*j* + 1), the term is again equally divided into the time derivatives of *P*_{ij}, *P*_{i+1,j}, *P*_{i, j +1}, and *P*_{i+1,j+1}.

The convected portion is mixed with all other grids, while the rest is assumed to be mixed with only grids in their vicinity; that is, the mixing ranges of *T* and *S* are in *T*_{i} − *T* ≤ *T* ≤ *T*_{i} + *T* and *S*_{j} − *S* ≤ *S* ≤ *S*_{j} + *S*. This choice reflects the condition that the properties vary gradually in a physical space except for the convected portion. Here the mixing for producing *P*_{ij} is represented by the mixing between a pair of the two grids symmetric about the (*i, j*) grid on the *T*–*S* plane (Fig. 6),

It is noted that this contribution is summed up for *P*_{ij} over the ranges of 0 ≤ *m* ≤ *T*/Δ*T* and 0 ≤ *n* ≤ *S*/Δ*S.*

The mixing coefficients *M*_{o} and *M*_{1} are estimated from scale analysis. The box is assumed to be filled with many blobs of water, each of which has a horizontal scale of *L.* Two blobs exchange water between them at a relative velocity of *u.* Hence, the water in one blob is replaced with the other blobs at a rate of *u*/*L.* The typical values of *L* = 100 km and *u* = 0.02 m s^{−1} give *M*_{o} = *M*_{1} = 2 × 10^{−7} s^{−1} as the reference values. Although the mixing range depends on actual mixing processes, the range is set to be (*T, S*) = (0.04, 0.01). Since the grid size on the *T*–*S* plane is chosen as (Δ*T,* Δ*S*) = (0.04, 0.01), the mixing in the nonconvected portion occurs between eight neighboring grids symmetric about a grid. Sensitivity to *M*_{o} and *M*_{1} will be examined in a parameter study.

## 4. Probability model results

The probability box model is used for examining the effects of heterogeneity on the equilibrium states that exist in the simple box model. A focus is given to the existence of the two equilibrium states, ES-2 (nonconvected) and ES-3 (convected); that is, whether both exist or only one of them exists. The probability model is first verified against the simple box model, where all parameters for the heterogeneity generation and the mixing are set zero. The lower box has (*T*_{o}, *S*_{o}) = (0, 0.4), and the initial value of the PDF is concentrated on one grid at (*T, S*) = (0, 0.1). The evolution of the spatially averaged state, ([*T*], [*S*]), is shown in Fig. 7 for 10 years along with the PDF. The probability model produces the averaged state similar to the simple model solution, converging to ES-2. The distribution indicates minor numerical effects, which would make *P* spread.

A parameter study is then carried out to examine the heterogeneity generation and the mixing with different values of *D*_{T}, *D*_{S}, *M*_{o}, and *M*_{1}. Once the mixing is included with *M*_{o} = *M*_{1} = 10^{−6} s^{−1} (five times the reference values), the PDF becomes more concentrated on one grid, as shown in Fig. 7. The heterogeneity generation is next introduced: *D*_{T} = 2 × 10^{−11} K^{2} s^{−1} and *D*_{S} = 10^{−12} ppt^{2} s^{−1} (20% of the reference values), and the mixing is set back to zero: *M*_{o} = *M*_{1} = 0. The averaged state shown in Fig. 7 is close to that in the simple box model, and the PDF spreads with the peak value greatly reduced.

The verification study has proved that the probability model reproduces the simple model solution, when the heterogeneity generation and the mixing are set zero. The parameter study is carried out with the mixing much stronger and the heterogeneity generation much weaker than the reference values. Their effects are significantly shown in the solutions. Hence, it is suggested that the PDF structure easily spreads due to the realistic heterogeneity generation, while the very intense mixing is required to converge the PDF effectively.

The existence of the two equilibrium states is examined with the lower-box property taken between (*T*_{o}, *S*_{o}) = (0, 0.3) and (0, 0.4). In the simple box model, only the ES-3 exists for *S*_{o} ≤ 0.35, whereas both the ES-2 and ES-3 can exist for *S*_{o} ≥ 0.36. The equilibrium values of the salinity are shown in Fig. 8. The salinity is in the vicinity of *S*_{o} for the ES-3, while it is smaller than *S*_{o} for the ES-2. As *S*_{o} increases, the ES-2 salinity would approach the value that is produced without vertical mixing.

The probability model is run for 30 years, in which a solution converges to an equilibrium state. For the case with the reference values of the heterogeneity generation and the mixing (*M*_{o} = *M*_{1} = 2 × 10^{−7} s^{−1}, *D*_{T} = 10^{−10} K^{2} s^{−1}, and *D*_{S} = 0.5 × 10^{−11} ppt^{2} s^{−1}), the equilibrium salinity and the probability of the convected portion, *P*_{IJ}, are shown in Fig. 8. Only one equilibrium state exists with any value of *S*_{o}, where *P*_{IJ} sharply varies in the range from *S*_{o} = 0.34 to 0.37. An additional calculation is carried out with *M*_{o} and *M*_{1} increased by a factor of 5. As shown by *P*_{IJ} in Fig. 8, two equilibrium states exist only around *S*_{o} = 0.37, where a large portion is convected at one state, and a small portion is convected at the other state. A major difference from the simple model is that the mostly convected state disappears, once *S*_{o} exceeds 0.38.

The evolution for the reference case with *S*_{o} = 0.37 is presented in Fig. 9 for 20 yr, along with the PDF at *t* = 20 yr. Here, the initial condition is given by the PDF concentrated on one grid either (0, 0.3) or (0, 0.37). It is of interest that the trajectories of ([*T*], [*S*]) are along that of the evolution in the simple model (Figs. 3 and 4). The evolutions of [*T*], [*S*], and *P*_{IJ} are also plotted in Fig. 10, showing convergence toward one equilibrium state in almost 20 yr.

The PDF distribution is compared with the *T*–*S* property observed by Johannessen et al. (1991). The ranges are ∼1°C and ∼0.2 ppt in the middle of the Greenland Sea. Thus, the observed distribution is wider along the *T* axis than the solution, while they are comparable along the *S* axis. This similarity is satisfactory in the present study, whose objective is to examine the effects of heterogeneity on the convected ocean using the conceptual model.

## 5. Discussion

The work in this paper is summarized as follows: the new box model, in which temperature and salinity are expressed by the probability distribution function (PDF), has been developed to examine convection in the subpolar ocean. The probability box model is built on the simple box model without PDF, which has only one active box for the upper layer of the subpolar ocean, receiving atmospheric forcing and interacting with the warmer box and the lower box. The simple box model potentially possesses a trimodal solution: the salinity-driven state (ES-1), the thermal-driven nonconvected state (ES-2), and the thermal-driven convected state (ES-3). The ES-2 is collapsed with the ES-3, as the lower-box density becomes close to the upper-box density.

The PDF is capable of describing horizontal heterogeneity in the upper layer. However, only a ratio of each water property is retained, while the information on its location is dropped. The mechanisms of heterogeneity generation and mixing are also parameterized. In the probability model with the heterogeneity generation and the mixing estimated for the realistic conditions, both of ES-2 and ES-3 do not appear; that is, either ES-2 or ES-3 exists. As the lower-box density varies, the transition between the ES-2 and the ES-3 occurs continuously.

Although modelers have been making efforts to include mesoscale variability and convection explicitly into a climate ocean model, available computer resources have not allowed them to do it. The model proposed in this paper has a highly simplified expression of spatial heterogeneity, whereas it may be useful for examining roles of mesoscale variability and convection in thermohaline circulation and climate change. The probability model has potential to be embedded in a low-resolution general circulation model (GCM) and a zonally averaged ocean model.

In a simple box model and a low-resolution GCM, the additional freshwater that is given to the models spreads over the sea surface and suppresses convection (e.g., Bryan 1986). Hence, the models are sensitive to a freshwater input. If mesoscale variability is modeled appropriately, there may be a weakly stratified portion favorable for convection. The convected portion acts as a favorable area for continuous convection. Thus, the system becomes more robust against the freshwater input. It is suggested that mesoscale variability should be treated carefully: the mesoscale variability is traditionally treated only as a mixing agent, whereas its heterogeneity itself should be included in a model.

The implication of the present model for the GIN Seas is discussed. Cessation of the deep-water formation in the GIN Seas has been observed in the last decade after the GSA, and may gradually modify the GSDW to higher temperature and salinity, and the lower density as well. Hence, the DSOW is reduced, and the weakened contribution to the NADW tends to increase its temperature and salinity. The simple box model predicts that the convection will be kept shut off as long as no heavy water is poured into the upper layer of the GIN Seas. The prediction by the probability model is different; that is, the convection may gradually recover in another decade. However, the ratio of the converted portion will be different, if the external forcing condition changed since the GSA.

## Acknowledgments

The work in this paper was started when the author worked at the Department of Fisheries and Oceans, Canada, with a partial funding support by the Canadian Federal Panel on Energy Research and Development. The author was supported also by Japanese Ministry of Education, Science, Sports and Culture. The author thanks J. Su, M. Yaremchuk, C. Griffiths, and W. Holland for fruitful discussion, and K. Ikeda for editorial reading.

## REFERENCES

### APPENDIX

#### Probability Distribution Function Arithmetic

The time derivatives of the spatial mean are derived from (3.1),

Equations (3.3) and the relationship (*T*_{i}, *S*_{j}) ≡ (*i*Δ*T*, *j*Δ*S*) lead to

Substitution of (3.2) for (*dT*/*dt*)_{ij} and (*dS*/*dt*)_{ij}, along with (3.1), yields (2.1), where *T* and *S* are replaced with [*T*] and [*S*].

The heterogeneity generation formulated by (3.4) is summed up for *P _{ij}*

This form in equivalent to the equation

for the PDF continuous on the *T*–*S* plane.

## Footnotes

* Current affiliation: Graduate School of Environmental Earth Science, Hokkaido University, Sapporo, Japan.

*Corresponding author address:* Dr. Motoyoshi Ikeda, Division of Ocean and Atmospheric Sciences, Graduate School of Environmental Earth Science, Hokkaido University, Sapporo 060, Japan.

Email: mikeda@eoas.hokudai.ac.jp