Insights on the non-linear solution of Munk’s ocean circulation theory from a rotating tank experiment Ocean and Coastal Research

At age 101, Walter Munk passed away in 2019. His groundbreaking discoveries will still guide and amaze oceanography students for years to come. Here, we perceive patterns in rotating tank with Munk’s circulation theory aided by a Lagrangian particles tracking algorithm and numerical modeling. From information captured by video, we track the trajectories of drifters, and then objectively mapped the streamfunction to obtain the mean circulation pattern. We were able to reproduce the wind-forced anticyclonic and asymmetric gyre, including the western boundary intensification and its retroflection. The latter phenomenon was predicted by the non-linear version of Munk’s model and observed in real subtropical gyres as small recirculation regions. We have configured two numerical model simulations mimicking the physical experiment, with linear and non-linear terms. The comparison between the numerical and physical experiments confirmed the effect of non-linear distortion of the gyre. Geophysical fluid dynamics is often hard to visualize, and counter-intuitive in a rotating system. We present this set of experiments as a tool for oceanography teaching. Besides studying general ocean circulation theories and observations through practical examples, this experiment provides an opportunity to develop basic image processing and geophysical modelling skills. ocean model.


INtrODUctION
The primary purpose of this article is to present a didactic fluid dynamics experiment in a rotating tank based on Marshall and Plumb's (2013) book to obtain not only qualitative evidence but a quantitative set of measurements to be further explored by graduate students. This laboratory experiment was designed to provide a simplified physical model of the mid-latitude, large-scale, wind-driven ocean circulation that retains the essential characteristics of its geophysical counterpart. As in most cases, (e.g.: Beesley et al., 2008), the focus is on the analogy between what is visually observed in the tank and how to link such patterns to the theoretical framework of simplified ocean models.
A lab demonstration by itself would allow hours of rich discussions and speculations about the nature of the movements involved. We go beyond this qualitative analysis and use the data recorded during the experiment to challenge students to extract quantitative information to substantiate their conclusions, regarding the similarity between the observed flow patterns on both tank and the ocean. A step further towards the refinement of the analysis includes the use of a general circulation numerical model.
Proficiency in the utilization of numerical models is an important and useful skill for aspiring academics for several reasons, especially due to lack of in situ data. Oceanic models present a wide range of applications, from idealized models built to isolate the effects of a single process at a time for hypothesis testing, to complex hindcast simulations employed in ocean state predictions. Here, we use the Regional Oceanic Modeling System (ROMS), a state of the art primitive equation solver oceanic model based on the Boussinesq approximation and the hydrostatic balance. First developed by Shchepetkin and McWilliams (2005), it has been widely used for scientific research.
The proposed experiment can be simplified for demonstrations to undergraduate students in the science, technology, engineering and mathematics areas. It can become more attractive to a more specialized audience of senior undergraduates or graduate students in physical oceanography and meteorology programs, or even to those in physics and engineering programs interested in experimental and computational fluid dynamics. Inspired by Marshall and Plumb (2013), we have increased the complexity of the experiment, extending (i) the theory that supports ocean and rotating tank observations, (ii) including quantitative image processing, (iii) data analysis techniques, and (iv) comparison with numerical experiments. As we further explore the quantitative aspects of the problem, we motivate students to develop computational abilities and apply them in problems that are useful for a wide range of oceanographic applications including satellite oceanography, Lagrangian trajectory analysis and regional modelling.

1.1.bAckgrOUND
We begin by presenting some key points relevant to the large-scale ocean circulation. These are intended to clarify which aspects of the theoretical framework are kept and which ones are lost due to the simplifications imposed to the physical model.

1.1.1.Subtropical gyres
Every major ocean basin has a subtropical gyre in which the currents, averaged over several years, show a smooth anticyclonic circulation. That circulation is mechanically forced by the wind stress against the sea surface in combination with the Coriolis force. In these subtropical regions, the wind gradually changes its direction from westward near 20º (also known as the Easterlies) to eastward near 40º of latitude (Westerlies), independent of the hemisphere. Due to the Coriolis effect, the wind forcing piles up water at the center of the gyre and introduces linear and angular momenta into the ocean. For simplicity, adiabatic conditions are assumed. The wind-induced circulation is counterclockwise in the Southern Hemisphere and clockwise in the Northern Hemisphere. In both cases, the presence of a bulge in the sea surface height corresponds to a region of atmospheric high-pressure in the mid-portion of the gyre (Figure 1).

Figure 1.
Colored contours of satellite-derived mean dynamic topography data from AVISO (Rio et al., 2018) where "H" indicates the location of atmospheric high-pressure centers. Black arrows represent an analytical function of the zonally and temporally averaged zonal wind stress that form the patterns historically named Easterlies (or Trades) and Westerlies.
Despite variations in the shape and size of the ocean basins due to the distribution of continental masses and bottom topography, there are three common aspects in mid-latitude circulation systems that are of dynamic origin. One is that western boundary currents (WBCs) are significantly faster than their eastern counterparts. This difference is apparent in Figure 1 for the WBC in the North Atlantic, where the Gulf Stream flowing along the US coast is much faster than the Azores Current on African/European side. Likewise in the South Atlantic, the Brazil Current is faster than the Benguela Current on the African side. A second point is that near the latitude of 35º, the poleward WBCs turn eastward quite abruptly and form what it is known as the WBC extensions. Finally, that change in direction is accompanied by a poleward overshoot that introduces perturbations into the flow. These perturbations are in the form of meanders and vortices that comprise a significant part of the energy of the system. These characteristics of subtropical gyres are a consequence of the conservation of angular momentum, or vorticity as oceanographers prefer. The steady-state balance can be summarized as follows: wind forcing accumulates mass in the center of the gyre via Ekman dynamic hydrostatics impose a horizontal pressure gradient force directed from the center to the borders of the gyre; the Coriolis force balances the pressure gradient and deflects the motion causing a (counter) clockwise basin-scale circulation in the (Southern) Northern Hemisphere. They were distilled in three classical ocean models from the 1950s: the Sverdrup, Stommel, and Munk models, that are discussed at length in Kundu and Cohen (2002), Vallis (2017) and other textbooks. The three models assume homogeneous density, hydrostatic balance, and a linear 2-D steady flow. Next, we highlight the most critical dynamical concepts with emphasis on the simplest model proposed by Sverdrup.

theOry
Physical processes can produce spatial and temporal phenomena that are distributed throughout a range of scales in the ocean, from millimeters and seconds to thousands of kilometers and decades, some with a time scale up to millennia. These movements are mathematically and physically described by the Navier-Stokes (Equation 1), an expression of momentum balance for fluids, and the equation of the conservation of mass. A detailed derivation is deferred to other studies, (e.g. : Pedlosky, 1986;Marshall and Plumb, 2013;Kundu and Cohen, 2002;Vallis, 2017). For the sake of simplicity, we choose the temporal scales adequate to our problem to reduce the momentum equation to: (1) where is the tridimensional velocity vector, is the gradient vector in Cartesian coordinates, Ω is the Earth's angular velocity, ρ is the mean density, p is the pressure, g is the gravity acceleration, and ν is the (turbulent) kinematic viscosity coefficient which is assumed constant.
A brief scale analysis, a concept accessible to undergraduate students, helps to evaluate the dominant balances as shown below: Asterisks denote non-dimensional variables with order 1 values. The uppercase letters represent typical scales of their respective lowercase dimensional variables. The horizontal scale L is used for both x and y, and the scales H and H v are used for z, the first being the depth of the gyre and the second the thickness of the viscous or Ekman layer. The vertical momentum component is assumed to be well represented by the hydrostatic balance ( ) and it is momentarily left out. To obtain a set of the most commonly defined non-dimensional numbers, the hydrostatic balance was used to scale the horizontal pressure gradient as ∆p = ρgH, make At this point we can build an a priori scaling argument to compare the physics of the tank and ocean. To that end, the scales listed in the upper half of Table 1 were used to compute the non-dimensional numbers listed on the lower half. Although not present in Equation 3, the Reynolds number, Re UH o = o , compares the acceleration and viscous terms, and it is often used as comparison between laminar and turbulent fluxes.
In Table 1, we chose L as the gyre width scale which is analogous to the diameter of the tank; H is the depth of the main thermocline to compare with that of the tank; H Re UH o = o in the ocean is the Ekman layer depth while in the tank it is the viscous layer thickness estimated from another experiment (Marshall and Plumb, 2013). The mean current velocity U can be roughly estimated from observations of the colored water motion in the tank; the value of the turbulent n 4 Polito et al. ν depends of the flux being analyzed thus we just stayed away from extreme values (Cushman-Roisin, 2011). The Rossby number, Ro, is a ratio between the acceleration and Coriolis terms and it is fairly small in both cases; the effects of rotation should not be overshadowed by transients in our experiment.
Although not as small as Ro, the Ekman number, Ek, that compares viscous and Coriolis terms is much less than one. That indicates that Ekman dynamics can be reproduced with the proposed experimental setup.
Bu can be interpreted in different ways (Cushman-Roisin, 2011;Vallis, 2017) and perhaps the most useful in this context is that it compares the length of the circulation with that traversed by the fastest gravity wave in one inertial period (i.e. the Rossby radius of deformation). Being order one means that any pressure perturbation will be influenced by rotation, and validates the application of a geostrophic framework to the case in focus. Finally, the values of Re indicate that although the ocean is turbulent, the flow in the interior of the tank is laminar, as the threshold Re is approximately 2000, (White and Corfield, 2006).
In Equation 4, we introduce the geostrophic velocity V g =(u g , v g ) that has only horizontal components, hence ∇ h =(∂ x , ∂ y ), and the vertical velocity (w) is zero. The Coriolis parameter f can be written as the scalar f=2Ω sin θ 0 , where θ 0 is a constant latitude. Notice that viscous layers do exist but they are very thin O[1] mm compared to the fluid depth in the tank O[10]cm and they are not explicitly resolved by geostrophy.
First let us consider the tank laying flat over the rotating table (θ 0 =90º). In this case, motion of a water column in any direction within the tank does not result in vertical stretching. Furthermore, the steady state or solid body rotation implies that the relative vorticity is constant over the whole tank, and in the non-inertial frame of reference it is null. This is an adequate representation of the so-called f−plane, f=2Ω, constant.
To mimic the models proposed by Sverdrup, Stommel, and Munk we need a key ingredient: the effect of the curvature of the planet. That cannot be neglected due to the meridional dimension of mid-latitude gyres. Physically it means that as a water parcel moves meridionally in the ocean it changes position relative to the axis of rotation of the planet. This changes the parcel's momentum of inertia and induces a change in angular velocity by conservation of angular momentum. We can mimic this effect by tilting the bottom of the rotating tank. This way, motion towards shallower or deeper sides of the tank implies in vertical stretching, that in turn changes the momentum of inertia of the water parcel. By the same fundamental principle, this will induce a change in angular velocity, similar to the planetary case. In this sense, the experiment with the tilted bottom rotating tank is analogous to the flat-bottom planetary beta-plane model. With that in mind, in the next sections we recall the geostrophic approximation on the β-plane to get to the simplified circulation models.

Geostrophy on the β-plane: the Sverdrup Balance
The f-plane becomes an inadequate approximation if the domain spans more than a few In this framework, our homogeneous slab of water can be safely assumed incompressible, thus we use the simplest version of the continuity equation ∇· u = 0. Because V g is horizontal, we can write: At this point it is crucial to notice that after taking the curl of the momentum equation we obtain an angular momentum equation, usually called a vorticity equation. When scaled by the layer depth, it becomes a potential vorticity equation. Equation 6 is a statement of conservation of potential vorticity in the absence of friction for a thin layer of homogeneous fluid over a sloping, rotating plane tangent to the surface of a spherical planet. The term on the lefthand side expresses the changes in potential vorticity as a particle is carried to the north or to the south. The term on the right is proportional to the vertical stretching or squeezing of the water column as it changes its momentum of inertia. At this scale, the compensation for vertical deformation is in the form of meridional motion over a sloping plane. Thus, in turn, we can simulate the β effect in laboratory by tilting the rotating tank giving it a linear topographic gradient.
In the ocean, the vertical stretching or squeezing of the water column can be introduced by the wind stress (τ) at surface or by the bottom friction. A balance is established between the Coriolis and the viscous terms, first and last of the right-hand side of Equation 1: In Equation 7, (u E , v E ) are the Ekman velocity components. From that balance applied at the surface or the bottom of the ocean, we could estimate the vertical velocity known as Ekman pumping (cf. Cushman-Roisin, 2011;Hendershott, 1987 and at the bottom: As Equations 8 and 9 came from the integration of Equation 7 (balance between the Coriolis and the viscous terms), it becomes clear that viscosity introduces vertical velocity. Viscosity is diffusion of momentum; in the ocean, it is effectively confined to thin layers of O[10 1 ]m. It has no direct means of transmission of momentum to the ocean interior layer, which is O[10 3 ]m thick. The mass convergence and divergence within that layer changes the pressure and generates a geostrophic circulation in the interior layer.
Next, we briefly discuss two more complex models that introduce boundary layers, bottom friction (Stommel, 1948), and lateral friction (Munk, 1950) in the Sverdrup balance.

Geostrophy on the β-plane-the Stommel and Munk models
Westward intensification was not explicit in Sverdrup's β-plane model, and to that effect Stommel (1948) introduced bottom friction through a linear dissipative term. This proved to be sufficient to reproduce the observed westward intensification of the wind driven currents. It is straightforward to show that geostrophic velocities (Equation 4) are non-divergent, and that allows the use of a scalar geostrophic streamfunction ψ to represent the flow. In terms of ψ, including the bottom friction as the last term, Equation 6 becomes: where , , and H >> h e .
Outside the Sverdrup interior, from Equation 10 we can depict the balance within a thin, meridional western boundary layer: , Stommel's solution for Equation 10 for a basin of width L connects the exponential decay at the border (x = 0) with a Sverdrup-like solution for the ocean interior (McWilliams, 2006): For a zonal, westward, sinusoidal wind stress of the form , with L y as the meridional length, the solution becomes: This solution represents the crowding of streamlines in the western boundaries of the subtropical gyres in the Stommel model.
A later addition to Equation 10 was the inclusion of the lateral friction by Munk (1950): where A h is the lateral viscous coefficient aforementioned. Munk (1950) assumed that the lateral friction of the form is much greater than bottom friction , thus neglecting the latter.
Munk's solution, including eastern and western boundaries (Vallis, 2017;Hendershott, 1987) is: Equation 15 reproduces the western boundary intensification as a meridional jet with a width of

ObjectIves
The practical exercise developed here provides an opportunity to elaborate on the main concepts of large-scale ocean circulation, in addition -and this is the novelty-to develop two important skills for physical oceanographers: data analysis and numerical modeling.
The material goal is to set up an experiment that reproduces the main aspects of the subtropical large scale ocean circulation, namely: • Wind Forcing: angular and linear momenta must be introduced at the upper boundary.  (Vallis, 2017). The main objective of this experiment is to substantiate laboratory observations with quantitative evidence. For that purpose, we chose a Lagrangian approach. However, a statistically significant number of particle trajectories should be obtained to reduce the noise introduced by perturbations caused from various sources such as from irregular blow patterns from fans, growth of small capillary-gravity waves from stress, convective patterns from heat exchange and evaporation, among others. These trajectories, or rather a statistically refined version of them, should allow us to address each one of the subtropical gyre circulation aspects mentioned above. The second point is a bit more subtle. If the flow that develops in the tank is indeed caused by the factors discussed in Section 1.2, a numerical model that solves these equations should provide similar results. Otherwise, some other physics outside the scope of the model may be important, e.g.: surface tension or heat exchange.

MethODs
Besides the rotating table with an adjustable speed and level, there are some materials necessary to reproduce the gyre circulation: • Acrylic tank (ours was a cylindrical tank with a 44 cm diameter but any tank shape will work); • Bubble level and a hose; • Wood planks; • 2 fans (computer coolers) and clamps to hold them in the tank; • Camera for still pictures or short videos; • Floats (Eppendorf tubes or paper confetti) and food coloring.

experIMeNtAl ActIvItIes
Our experiment is divided into three parts: laboratory procedure in the rotating tank, the particle trajectory analysis, and numerical modeling.

The rotating tank preparation
The rotating table is leveled and the camera is fixed to a support arm just above the tank. The camera records the images in a non-inertial reference frame because it rotates together with the tank. The tank is centered on the table and it is tilted to about 7º using wood planks to simulate a meridional gradient of the planetary vorticity (see Figure 3 and the discussion surrounding Equations 4 and 5). The tank is filled with water up to a level sufficient to fit the fans at the surface. A North sign fixed in the tank would be convenient for interpreting the circulation patterns later in the pictures. This would be a good opportunity to encourage the students to discuss the physical meaning of tilting the tank, and what position in the tank represents the North. The linear variation of depth is used here to simulate the β-effect in geophysical fluids. We see from Equation 6 that the topographic β is inversely proportional to the changes in the planetary vorticity. Thus, the North in this experiment coincides with the shallower side of the tank. Next, the fans are attached on opposite sides of the tank with clamps to produce a wind pattern that can transfer negative vorticity movement to the water. It is desirable to have fans adjusted to blow approximately parallel to the surface. At this point, the students are invited to discuss the correct position and facing direction of the fans to generate the desired wind circulation.

The Laboratory Experiment
The experimental procedure is summarized in Figure 4, upper half. With the fans off, the rotating table is set to rotate counterclockwise at approximately 10 rpm until it reaches the rigid body rotation, after approximately 25 minutes. At this stage, confetti put at surface should not be moving in the rotating frame of reference. The fans and the camera can be turned on. The camera should be previously set to record a short movie (our approach for this experiment), or to take pictures at predetermined time intervals (for example, every 5 s). To identify the main current pathways, we used water-based food coloring. Contrasting dye can be carefully introduced at the surface, for instance in the NE and SW sectors of the tank. Once some patterns are visible, floats are deployed to work as targets to have their velocity estimated by the particle tracking algorithms. As a suggestion, stable floats can be made with 0.2 ml Eppendorf tubes partially filled with water to balance their buoyancy and black confetti glued on the top cap (to facilitate tracking). The goal of this part is to obtain a collection of pictures or a movie.

pArtIcle trAjectOry ANAlysIs
Quantitative data on the float movements are obtained using particle detection and trajectory estimating algorithms applied to the images recorded by the camera. The steps for data processing are summarized in Figure 4, lower half.
Either working with a series of still pictures or a movie, the first step is to extract the frames and convert them into some popular image format, e.g. PNG, JPEG, or TIFF. That can be done by using an open source software such as ffmpeg 1 . A frame rate of 0.5 frames per second for the conversion of movie frames works fine; as for the case of shooting of still pictures, an interval of 5 seconds is sufficient. After the images have been converted to grayscale, a particle tracking algorithm should be used to extract the data. There are some open source software packages available in public repositories. We used Trackpy 2 software (Allan et al., 2016) for particle detection and trajectory determination. This software requires as input the size (in pixels, hereafter px) of the observed particles, e.g.: the radius of the confetti.
For each identified feature, the software produces a parameter called mass that quantifies its luminosity density peak. A threshold mass value can be defined to filter out specks that qualify as spurious data. That part of the procedure partially removes undesirable data resulting from light reflection, the shadows of the clamps and fans, etc. The estimation of the trajectory of an individual float is based on the path length (in px) traversed from one frame to the next. We set a value of 200 px as the displacement limit between frames. Based on the observation of trajectories, sometimes it is necessary to set a maximum value of frames for a particle that vanished to reappear following the same path. We used a value of 10 frames.
Even after all image processing procedures applied up to the steps depicted in the last row of Figure 4, we observed that some of the wrongly identified trajectories represented almost inert objects. These are due to the presence of shadows, refraction and reflection effects, such as the bright refraction areas off the tank's edge and shadows of the fans. These light artifacts do not move much in the rotating frame of reference where the camera is mounted and their geometry does not resemble that of our Lagrangian floats. Therefore, quasi-stationary objects moving less than 2 px between frames were removed. There were also trajectories that developed during a short period of time and disappeared. These are not associated to any of the floats but again to light reflecting off surface ripples, thus in the last step we consider only trajectories that lasted longer than 5 frames. Finally, the velocity from the selected trajectories was estimated by dividing the displacement between two consecutive frames by its time lapse. Since we are interested in obtaining the large-scale circulation pattern, after converting the velocity values to ms -1 , we bin averaged the data to smooth the spatial (60 px) and temporal (6 s) scales (i.e., the inertial period, since our tank rotates at 10 rpm).

2.2.1.Objective Analysis
This experiment presents a great opportunity to deal with analysis of irregularly-distributed noisy data. Shear instability on the western boundary shed eddies that could influence the whole domain. Therefore, to analyze each identified trajectory individually in a noisy velocity field could be misleading. We minimized those effects by mapping the velocity and streamfunction from the Lagrangian analysis of the floats' trajectories. For that, we performed an Objective Analysis (OA), an optimum interpolation technique first used in oceanographic data by Bretherton et al. (1976) and later improved by Carter and Robinson (1987). A gridded field can be produced from a linear combination of the observed data (φ i ), weighted by a correlation function (C ij ): Equation 16 was obtained following Carter and Robinson (1987)'s methodology, which assumes that it fits an isotropic Gaussian function, such that: where 2 e is the variance of the random error, is the distance between the observation points, and l c is the radial correlation length scale.
To determine the form of the correlation function, all observations within 35 px from each other were grouped together and the correlation between them was 1 Video processing software under LGPL licensing https://www.ffmpeg.org 2 Particle tracking algorithm written in Python http://github.com/soft-matter/trackpy Polito et al. calculated. This procedure was repeated progressively by increasing the lag between observations. The correlation distribution ( Figure 5A) allowed us to determine the l c ≃ 5 cm (94.240 px) and the 2 e ≃ 0.116 used in Equation

17
. Figure 5B shows the individual particle's trajectory. Using the correlation function properties, we mapped the OA error in a gridded field ( Figure 5C), and used it to mask out the regions where the interpolation error might affect the velocity results. We applied the noslip (Dirichlet) boundary condition at the limits of the domain in Figure 5C. In practice, a circular mask was applied to the objective analysis domain to mimic the tank limits, enforcing zero velocity at and beyond the tank boundary.

The Numerical Experiment
We ran the ROMS-RUTGERS ocean circulation model because it is a popular, up-to-date, and easyaccess ROMS version (Shchepetkin and McWilliams, 2005). The model comes in three versions (RUTGERS, AGRIF, UCLA) with the RUTGERS and AGRIF being open source. Although the AGRIF version is the easiest one to use, we have chosen RUTGERS version because they maintain an active discussion community and perform frequent updates to its source code. We created a 44 cm square grid with an outer circular mask. The spatial resolution of the grid was 2.2 × 10 -1 cm, with a time step of 10 −2 s. These values are set to avoid violating the Courant-Friedrichs-Levich condition, necessary for the numerical convergence of the model. The Coriolis parameter of 2.1 s -1 was calculated from the rotation rate. The lateral boundaries of the domain are closed so that the only forcing comes from the wind at the surface, set to increase radially from the center of the grid. The bathymetry varies linearly from 21.1 cm to 15.7 cm on the meridional axis, with the deepest point on the southern portion of the grid, approximating what is shown in Figure 3. We employed a logarithmic drag and set the water temperature to 24ºC in the simulation. We let the simulation run for one hour, and saved an output at each half minute, generating a total of 120 snapshots. Further details of the model configuration are deferred to the supplementary material.

resUlts AND DIscUssION
The OA was calibrated with the autocorrelation function and the associated error of the velocity field. The resultant interpolated field is shown in Figure 6A, where masked areas correspond to formal errors greater than 0.1 in Figure 5C. The east-west asymmetry of the wind-forced gyre is conspicuous. The WBC flows northward, and after crossing the halfway point in the tank, the current retroflects and forms a large anticyclonic recirculation gyre. The OA yields an error-constrained smooth streamfunction of the mean circulation, where the northward flowing WBC deviates to the east and meanders as it separates from the western boundary. Taking the radius of the recirculation gyre as 5 cm (same as the estimated correlation length scale) and a typical velocity scale in the western boundary (1 cm s −1 ), the Rossby number ( ) f g is O[10 -1 ], which indicates that the nonlinear terms are small but cannot be neglected for this feature. Indeed, the recirculation pattern clearly seen in the streamfunction field resembles the non-linear/ no-slip Munk solution for the subtropical large-scale circulation, presented by Vallis (2017), and adapted here in Figure 6B. As explained in the introduction, the goal of this experiment was to represent the same features described by the linear theory, but it turns out that the results indicate that the non-linear terms are important to describe the overshooting lobe observed in the NW part of the rotating tank gyre. This is also an opportunity to discuss with the students the role of the non-linearity on geophysical fluids.
The two approaches employed in our data analysis, the Lagrangian trajectories ( Figure 5B) and the OA streamfunction ( Figure 6A), represent the observed circulation patterns. The insight that students gain undoubtedly goes beyond the main objective of this activity that is to connect the observations with the theory. They get hands-on experience with (i) a particle-tracking software based on image analysis and (ii) a reasonably sophisticated statistical interpolation method. The programming and data analysis skills acquired in this exercise are useful for a broad range of geophysical applications.
Qualitatively, the numerical model simulations had all the essential features observed in the rotating table experiment: forcing from the upper boundary, the Northern Hemisphere rotation and β-effect, and as expected, an asymmetric gyre produced an intense WBC that retroflects at a given latitude and non-linear distortion in the form of a big meander downstream of the retroflection. Figure 7A shows the results from the linear experiment. The positive ψ values on the western boundary indicate high pressure on that side of the tank and consequently a clockwise circulation (anticyclonic for the Northern Hemisphere), in agreement with the "Stommel solution" (Stommel, 1948), presented in Equations 10 to 13. The WBC intensification is also inferred by the convergence of streamlines on the western portion of the simulated tank. The linear experiment output is an "order zero" representation of the most basic results found in our physical experiment. Despite the difference in the geometry of the domain, there is a striking similarity between our OA solution in Figure 6A and the "Munk solution" ( Figure 6B; see also Equation 14 and Vallis, 2017). We performed a non-linear numerical simulation and the outputs are presented on Figure 7B. The addition of a non-linear component in the modelwhich has to be accounted for when working in a small tank such as ours-greatly increased the similarity between the OA observations ( Figure 6A) and the non-linear run of the ROMS model ( Figure  6B). The perturbations of the streamlines resemble the observations in the ocean and in the tank, stressing the importance of non-linear effects on the large-scale circulation.
Given the simplifications of Munk's model, a rigorous quantitative match between model and experiment is not expected. Differences are expected as consequence of (i) ignoring the latent heat exchange in the top boundary, (ii) improper modeling of the wind stress horizontal distribution and vertical momentum transfer at the surface, (iii) short duration of the experiment allows for transient eddies to influence the results, and (iv) other minor effects. With that in mind, the velocity magnitude can be inferred from both Figures 6A and 7B from the r T T} ratio, where r is the distance along a section perpendicular to } contours in the western boundary current region. In both cases one arrives at a value on the order of 0.01 m.s -1 . Although the location of the gyre center is slightly more to the North in the model, the width of the Munk layer L M is similar. This length can be estimated for the tank, with few considerations. From the theory in the tank, β is simulated with the bottom slope α, such that Using A h = 10 -6 m 2 s -1 results in a Munk layer of L M = 1 cm, similar to the value of the observed WBC in the tank ( Figure 5B) and the same order of magnitude of the WBC in the numerical model ( Figure 7B). The nonlinearity parameter S L U 2 b = mentioned in Figure 6B can be estimated for our experiment. Using the same b=1.6 (m.s) -1 , U=0.01 m.s -1 and L as the diameter of the tank, one arrives at S=0.13 (the case shown in Figure 6A), similar to the value used by Vallis (2017).

cONclUDINg reMArks
In this paper we proposed a multi-faceted didactic approach to study a major feature of the ocean circulation: the wind-driven subtropical gyre. A data-analysis tool is assembled to process and get useful information from a video recorded geophysical fluid dynamics experiment, in an Polito et al. analogous fashion to satellite remote sensing data processing. In the proposed activity, we reproduced the wind-forced, anticyclonic, asymmetric-clockwise gyre with the western boundary intensification, both numerically and in a rotating frame. We also observed the retroflection of the intensified western jet to conserve vorticity, as observed in the recirculation regions of the WBC of the subtropical gyres. The comparison between both the numerical and laboratory experiments confirmed the expected effect of non-linearity distortion within the gyre.
From an educational perspective, besides reproducing general ocean circulation theories and observations through practical examples, this experiment provides an opportunity to further understand non-linear effects on the large-scale ocean circulation.
Polito et al.

sUppleMeNtAry MAterIAl the NUMerIcAl MODel set-Up
This material supplements the manuscript "Insights of the non-linear solution of Munk's ocean circulation theory from a rotating tank experiment". Here we describe in more detail the Regional Oceanic Modeling System (ROMS) configuration necessary to reproduce the numerical experiments. There is extensive documentation in the model website and forum that explains in detail the installation procedure. We assume that the reader already has a working version of the ROMS-RUTGERS model installed.
There are two main ways to set up the model: (i) configuring in its own files in the User folder of the model or (ii) directly changing the model source code. We will do the latter, for the sake of simplicity. The first step is to define a project name for this run of the model. In the Makefile file, this name will be defined by the ROMS APPLICATION variable. This step is crucial for the rest of the configuration since it will be used to reference the initial condition and forcing in the analytical files.
The next step is the creation of a grid, and again, there are two ways of doing that: one is to change the file ana grid.h, and the other is to create a NetCDF grid file. We took the second approach, and for that we need to input four fields to be read by the model: grid coordinates, mask, bathymetry, and Coriolis frequency. We used Python to create our grid with the pygridgen library. We used a circular tank with 0.44 m in diameter, thus our numerical domain is 0.44 m wide in both x and y dimensions. We used a grid spacing of 0.022 m to guarantee a proper resolution to visualize the boundary current. We also included a masking field in our domain. In that field, we set all values outside an inscribed circular region to 0. In other words, supposing the grid was centered at 0, then all values where |x| ≥ 0.22 m were set as 0, where| • | is the euclidean norm, and x = (x, y) are the discrete grid coordinates in meters. For the bathymetric field we defined the following function: where L = 0.44 m, h max = 0.211 m, and h min = 0.16 m which was defined as being the minimum depth where the tank achieved a 7º angle tilt. With that definition we guarantee a linear decrease of the depth in the y axis, as was configure in the field experiment. Finally, the Coriolis frequency was set to 2.1 s −1 in all grid points as to simulate a 10 rpm rotation in the f-plane. All four fields (grid points, mask, bathymetry, Coriolis) were saved in a NetCDF file to be read by ROMS. The Octant and Pygridgen libraries are openly available to the public at GitHub. The following code was used for creating the grid: import numpy as np import pygridgen as pg import octant . roms as rms x = np. linspace (0, 0.44, 201) # x is dimension y = np. linspace (0, 0.44, 201 # y is dimension xg, yg = np. meshgrid (x , y ) # creating spatial grid grd = pg. CGrid (xg, yg) # creating roms grid in Arakawa−C # creating bathymetry h = − (0.211 − 0.089) / 0.44 * grd . y rho + 0.211 # defining radius r = np . s qrt ((grd . x rho − 0 . 22) ** 2 + (grd . y rho − 0 . 2 2) ** 2) grd . mask [ r > 0 .22] = 0 # creating mask # defining Coriolis parm. grd . f = 2 * (2 * np . pi /( 60.0 / 10)) * np . ones (grd . x rho.shape) grd . h = h # Input Bathymetry files = 'roms grid . nc' #pathtosavefile rms.write grd (grd, filename = files) #Savinggridinnetcdf The next step is to define which physical parameters should be used by the model (e.g.: longwave and shortwave radiation, temperature and salinity advection, momentum advec-tion, etc.). In ROMS terms, those are defined as flags, and they are defined in the Include folder. In that, the file cppdefs.h shows all available flags. We set the following: • To test the importance of nonlinearity in our experiment, we ran two models, one including the advection flag, the other excluding it. It should also be noted that, if the reader chooses to use the grid from the model itself rather than creating a separate one, the ana grid and ana mask flags should also be included, and their respective files modified.
Since we set the analytical surface forcing, we needed to change the respective file in the Functionals folder. For that, we modified the ana smflux.h to include the wind stress. Our laboratory experiment used two fans at the zonal edges of the tank, positioned lower than the border to confine most of the air motion. To simulate this, we defined the following wind stress: where W = 3.89 10 −5 N m −2 is maximum wind stress, located at the boundaries. This value was calculated using the equation: ( 3) where U=0.0057 m/s is the wind speed extracted from the manufacturer's information, ρ air = 1 kg/m 3 is the air density, and C d =1.2 is the drag coefficient extracted from (Large and Pond, 1981), we then arrived to our wind stress amplitude. Although we have no interest in the analysis of net heat or bottom temperature fluxes, these parameters should be included because the model requires them. From the analytical parameters, however, the only files that should be edited are the ANA INITIAL, and ANA SMFLUX. For the ANA INITIAL, it is written inside an if-else condition, as follows:
The last step is to configure the external file in the External folder to run the model. (Note that this file includes several parameters which we are not dwelling with.) For our model, we used the upwelling test case file as a template and changed these parameters: • The project name "MyAppCPP" has to be the same one used in the Makefile, and in the analytical surface forcing; • The grid dimension parameters; • The momentum boundary conditions, which in our experiment were set to "closed"; • The slip conditions (gamma2), which in our experiment were set to "no-slip", i.e.-1; • The time stepping, which has to be set as to not violate the Courant-Friedrichs-Lewy (CFL) condition; • The output parameters (e.g. variables to be saved, saving frequency); • The grid path to reference our grid file.
ROMS is able to resolve much more complicated problems. Because we are dealing with a highly simplified configuration of the model, there are a large number of parameters left with their default values. These are not important to our particular purpose, however knowing their purpose may be useful for future developments. This exercise is intended as a first step in the use of numerical modelling.