Treatment of rivers
Treatment of rivers
Is it our consensus that the best way to put freshwater in is as a surface plume instead of as a well-mixed layer?
I would have thought to put them in as a well-mixed layer with freshwater spread throughout all vertical levels.
Courtney and Aaron
I would have thought to put them in as a well-mixed layer with freshwater spread throughout all vertical levels.
Courtney and Aaron
I remember somebody (either John Warner or Hernan) saying the rivers should be put into the top layer. Maybe they can help us with a primer on how to do rivers.
Chris
csherwood@usgs.gov
Chris
csherwood@usgs.gov
Last edited by csherwood on Fri Apr 02, 2004 7:00 pm, edited 2 times in total.
It was my experience with ROMS 12 months ago that putting the rivers into the top layer worked best in the respect that it induced less over/under-shoot of the salinity solution. I was getting salinity values of 40 PSU in the neighbourhood of the river mouth! Since then, some bugs have been fixed in the code (associated with the application of the masks) and I have not rechecked this issue.
An argument has been made that it's also important that the river properties of salinity and temperature be prescribed, not inferred with what amounts to a gradient boundary condition (we allowed this option for temperature), because this leads to an unconditionally unstable advection step.
To me it's more rational to put the water in the top layer because distributing it uniformly in the vertical means than on every time step you are imposing the internal dam burst problem and the model has to try to adjust a vertical isopycnal to a horizontally stratified salt wedge by internal waves. Effectively, the hydrostatic approximation is violated on every step. This is a very hand waving argument, but if the grid cells are 8 km across (as they are in my application) then it seems to make sense to impose an inflow that looks something like how I expect the ocean to adjust physically over the 8 km.
I did the math something like this: If the cell is dx wide, h deep, and I put the flow in through the topmost fraction f of the water column, then the average inflow velocity is u = Q/(f*h*dx) m/s (for a river transport of Q). If your cells are small and shallow, then u could be unrealistically large for your river, so you might want to have f be a substantial fraction of the water column. In my case, making f=1/N for h=20 m still gave wicked small u, so I was unconcerned about that.
Enrique might have a counter argument to this for the case where the depth at the coast is quite large and an inflow in the top cell only forces a bottom detached plume, which might be undesirable.
This needs a systematic study in a small idealized test case. Any volunteers?
John
wilkin@imcs.rutgers.edu
An argument has been made that it's also important that the river properties of salinity and temperature be prescribed, not inferred with what amounts to a gradient boundary condition (we allowed this option for temperature), because this leads to an unconditionally unstable advection step.
To me it's more rational to put the water in the top layer because distributing it uniformly in the vertical means than on every time step you are imposing the internal dam burst problem and the model has to try to adjust a vertical isopycnal to a horizontally stratified salt wedge by internal waves. Effectively, the hydrostatic approximation is violated on every step. This is a very hand waving argument, but if the grid cells are 8 km across (as they are in my application) then it seems to make sense to impose an inflow that looks something like how I expect the ocean to adjust physically over the 8 km.
I did the math something like this: If the cell is dx wide, h deep, and I put the flow in through the topmost fraction f of the water column, then the average inflow velocity is u = Q/(f*h*dx) m/s (for a river transport of Q). If your cells are small and shallow, then u could be unrealistically large for your river, so you might want to have f be a substantial fraction of the water column. In my case, making f=1/N for h=20 m still gave wicked small u, so I was unconcerned about that.
Enrique might have a counter argument to this for the case where the depth at the coast is quite large and an inflow in the top cell only forces a bottom detached plume, which might be undesirable.
This needs a systematic study in a small idealized test case. Any volunteers?
John
wilkin@imcs.rutgers.edu
Just to expand on the counter argument...
In a situation like the North American West Coast, in particular the Columbia river, putting all the flowin the top layer causes a plume to extend (many) hundreds of kilometers in to the Pacific, capping the surface with a fresh layer that shuts down any mixing that needs to take place. Distributing the inflow over various levels has created a more realistic looking plume, though how to do the distribution is so far a bit of an art. Lacking a proper estuary model, we do need to better understand how the model is handling the river inputs. I am willing to put some time into this this summer since this issue has not gone away for some time now. It might help to make a list of the possible variables.
As John started saying the variables may be:
1. Horizontal (and vertical) resolution
2. Inflow rate
3. Depth at source
4. Tracers specification
5. Bottom drag?
Enrique
In a situation like the North American West Coast, in particular the Columbia river, putting all the flowin the top layer causes a plume to extend (many) hundreds of kilometers in to the Pacific, capping the surface with a fresh layer that shuts down any mixing that needs to take place. Distributing the inflow over various levels has created a more realistic looking plume, though how to do the distribution is so far a bit of an art. Lacking a proper estuary model, we do need to better understand how the model is handling the river inputs. I am willing to put some time into this this summer since this issue has not gone away for some time now. It might help to make a list of the possible variables.
As John started saying the variables may be:
1. Horizontal (and vertical) resolution
2. Inflow rate
3. Depth at source
4. Tracers specification
5. Bottom drag?
Enrique
Last edited by enrique on Fri Apr 02, 2004 7:01 pm, edited 1 time in total.
Rivers are a boundary condition within the model. As such, they should be formed just like any other boundary condition -- in a way that is numerically sound and physically consistent.
That said, a river boundary condition should look like river flow, especially if this is put at the head of some sort of estuary. You could have some sort of boundary layer (quadratic profile) or constant stress (linear profile) that is zero at the sea floor, and whatever it needs to be at the surface to maintain the mass flux. I typically use the linear profile.
There are times you might choose to only put the fresh water in some surface layers. For example, I put the Mississippi outflow in the upper 2 m of the water column (over about 5 points), because that is really what it looks like. Again, it is a linear profile with the flow being zero at and below 2 m.
You should never, ever put the flow in the upper layer alone. Numerically, you need to resolve the vertical structure of the outflow, meaning that the outflow should vary over at least a few vertical points. Otherwise, there is unrealistically high gradients (and shear) in this region, and numerical instabilities will form. This is also a completely unphysical representation of what a buoyant input would look like. Bad, bad, bad...
Rob
That said, a river boundary condition should look like river flow, especially if this is put at the head of some sort of estuary. You could have some sort of boundary layer (quadratic profile) or constant stress (linear profile) that is zero at the sea floor, and whatever it needs to be at the surface to maintain the mass flux. I typically use the linear profile.
There are times you might choose to only put the fresh water in some surface layers. For example, I put the Mississippi outflow in the upper 2 m of the water column (over about 5 points), because that is really what it looks like. Again, it is a linear profile with the flow being zero at and below 2 m.
You should never, ever put the flow in the upper layer alone. Numerically, you need to resolve the vertical structure of the outflow, meaning that the outflow should vary over at least a few vertical points. Otherwise, there is unrealistically high gradients (and shear) in this region, and numerical instabilities will form. This is also a completely unphysical representation of what a buoyant input would look like. Bad, bad, bad...
Rob
It's interesting to hear other user's experiences here, but I would caution against simply doing what someone else does without checking closely... hence my call for a methodical test.
The Mississippi is certainly an exceptional river in terms of volume transport, so Rob Hetland's comment that putting the inflow into as little as 2 m of the upper water column works well is illuminating.
Rob doesn't say, but if he's working with grid cells as large a 5 km, then my logic in the previous message is borne out. Lifting some nifty factoids from
http://www.nps.gov/miss/features/factoids/
and the google calculator (to convert pesky US units to metric) we see that a mean Mississippi transport of Q = 600,000 feet^3/s (1.7e4 m^3/s) distributed across a 5 km wide cell over 2 m in the vertical would be entering at, on average, 1.7e4/(2*5000) = 1.7 m/s. This seems pretty quick, especially if the profile is linear over the 2 m implying the surface current is 3.4 m/s, but the other factoid I found was that at New Orleans, on 2/24/2003, the speed of the river was 3 miles per hour which is 1.3 m/s (not too far from the 1.7 m/s average). So this works nicely for inputting a whopper of a river like the Mississippi.
But consider a more modest river like the Hudson at 600 m^s/3 on average. Entering a 5 km wide cell in a 2 m layer gives a trickle of 6 cm/s. If the model water depth were 5 m and there were 20 vertical levels, i.e. dz = 0.25 m, putting the flow into the top level only would imply u = 25 cm/s, which is not outrageous. So I don't agree with Rob's assertion that ...
John
wilkin@imcs.rutgers.edu
The Mississippi is certainly an exceptional river in terms of volume transport, so Rob Hetland's comment that putting the inflow into as little as 2 m of the upper water column works well is illuminating.
Rob doesn't say, but if he's working with grid cells as large a 5 km, then my logic in the previous message is borne out. Lifting some nifty factoids from
http://www.nps.gov/miss/features/factoids/
and the google calculator (to convert pesky US units to metric) we see that a mean Mississippi transport of Q = 600,000 feet^3/s (1.7e4 m^3/s) distributed across a 5 km wide cell over 2 m in the vertical would be entering at, on average, 1.7e4/(2*5000) = 1.7 m/s. This seems pretty quick, especially if the profile is linear over the 2 m implying the surface current is 3.4 m/s, but the other factoid I found was that at New Orleans, on 2/24/2003, the speed of the river was 3 miles per hour which is 1.3 m/s (not too far from the 1.7 m/s average). So this works nicely for inputting a whopper of a river like the Mississippi.
But consider a more modest river like the Hudson at 600 m^s/3 on average. Entering a 5 km wide cell in a 2 m layer gives a trickle of 6 cm/s. If the model water depth were 5 m and there were 20 vertical levels, i.e. dz = 0.25 m, putting the flow into the top level only would imply u = 25 cm/s, which is not outrageous. So I don't agree with Rob's assertion that ...
I stand by my heuristic argument that it is better to err on the side of entering the flow in a thin layer to avoid the endless internal dam burst problem, which is not only unphysical but formally violates the hydrostatic assumption that underpins the ROMS governing equations.You should never, ever put the flow in the upper layer alone. Numerically, you need to resolve the vertical structure ofthe outflow, meaning that the outflow should vary over at leasta few vertical points. Otherwise, there is unrealistically high gradients (and shear) in this region, and numerical instabilities will form. This is also a completely unphysical representation of what a buoyant input would look like. Bad, bad, bad...
John
wilkin@imcs.rutgers.edu
The Mississippi Delta is like a leaky sponge, with a few cuts in it. Southwest Pass (the major outflow) is resolved with three horizontal points, (slightly less than 1 km wide each), and it carries about 1/3 of the total flow. Maximum transports are about 4000 m3/s, and velocities in the vicinity of the plume can be greater than 2 m/s, with supercritical flow. I have 43 sources all around the delta to get the outflow looking right.
This is obviously a very different scenario than a model with 5 km resolution. In this case, the strong mixing and advection that occurs near the mouth of a narrow estuary is completely unresolved -- mixing is done completely by the wind (and to a much lesser extent, other large scale processes) and the flow is always subcritical.
I would propose that the outflow should be mixed somewhat already in these large-scale runs. It would be better to explicitly parameterize the sub-gridscale mixing by increasing the salinity (and, correspondingly, the transport) of the outflow. You can expect the model with 5 km resolution to get it even close to right, why not admit that this is a sub-gridscale process, and parameterize it in a way that makes sence. You might need to take out salt and mass in the lower layer if you are worried about salt conservation (and you can do this with a 'negative' source, but you can't yet make the reference salinity of the negative source be determined on the fly by the ambient water on the shelf. This, obviously, could get out of hand if there was a mismatch between what is there and what is sucked out..). I don't think there are any perfect, or even good, solutions at the moment; this is a wide open problem right now...
However, I stand my my assertion that nothing good can come of putting the fresh water in the upper layer alone. I am less concerned about unphysical numbers as I am about numerical instability when anything is resolved with only one grid cell. The first rule of numerics is that features need to be resolved. Even our fancy, high-order advection schemes need a few points to resolve a feature without terrible numerical dispersion. Why not put the outflow in the upper three cells, with the largest outflow at the surface? It seems to me that, physically, this is very similar, and numerically much more sound.
Rob
This is obviously a very different scenario than a model with 5 km resolution. In this case, the strong mixing and advection that occurs near the mouth of a narrow estuary is completely unresolved -- mixing is done completely by the wind (and to a much lesser extent, other large scale processes) and the flow is always subcritical.
I would propose that the outflow should be mixed somewhat already in these large-scale runs. It would be better to explicitly parameterize the sub-gridscale mixing by increasing the salinity (and, correspondingly, the transport) of the outflow. You can expect the model with 5 km resolution to get it even close to right, why not admit that this is a sub-gridscale process, and parameterize it in a way that makes sence. You might need to take out salt and mass in the lower layer if you are worried about salt conservation (and you can do this with a 'negative' source, but you can't yet make the reference salinity of the negative source be determined on the fly by the ambient water on the shelf. This, obviously, could get out of hand if there was a mismatch between what is there and what is sucked out..). I don't think there are any perfect, or even good, solutions at the moment; this is a wide open problem right now...
However, I stand my my assertion that nothing good can come of putting the fresh water in the upper layer alone. I am less concerned about unphysical numbers as I am about numerical instability when anything is resolved with only one grid cell. The first rule of numerics is that features need to be resolved. Even our fancy, high-order advection schemes need a few points to resolve a feature without terrible numerical dispersion. Why not put the outflow in the upper three cells, with the largest outflow at the surface? It seems to me that, physically, this is very similar, and numerically much more sound.
Rob
About not putting all the flow into the top layer and the model numerics:
We may have some evidence for not putting it all in the surface layer. When the river discharges are bumped up unexplainable (at least by me) floating point and overflow errors seem to miraculously occur, sometimes. The only explanation I could come up with is that the model does not like that high of a gradient between all the fresh and salt water, which goes along with :
"I am less concerned about unphysical numbers as I am about numerical instability when anything is resolved with only one grid cell. The first rule of numerics is that features need to be resolved. Even our fancy, high-order advection schemes need a few points to resolve a feature without terrible numerical dispersion."
Aaron
We may have some evidence for not putting it all in the surface layer. When the river discharges are bumped up unexplainable (at least by me) floating point and overflow errors seem to miraculously occur, sometimes. The only explanation I could come up with is that the model does not like that high of a gradient between all the fresh and salt water, which goes along with :
"I am less concerned about unphysical numbers as I am about numerical instability when anything is resolved with only one grid cell. The first rule of numerics is that features need to be resolved. Even our fancy, high-order advection schemes need a few points to resolve a feature without terrible numerical dispersion."
Aaron
John Wilkin and I had a discussion this week about a boundary condition for the Hudson's outflow into a model of the New York Bight. One issue that came up is that the salt balance at the tidal time scale in many (most?) estuaries is not in steady state. This means that the flux of fresh water out of the system will modulate over the tidal month-- even if the freshwater discharge into the estuary is held constant.
Is this important to the plume? Don't know-- but since the Spring/Neap cycle can (in some estuaries) dramatically modulate both the volume transport and salinity of the outflow, it must have at least some near-field effect to the coastal current.
How far down stream does this impact the coastal current? don't know that either-- but idealized runs with an estuary undergoing spring neap cycles should be assess to answer this. In the real world, spring neap modulation at the plume's source is likely to be eventually obscured by winds.
Bob Chant
Is this important to the plume? Don't know-- but since the Spring/Neap cycle can (in some estuaries) dramatically modulate both the volume transport and salinity of the outflow, it must have at least some near-field effect to the coastal current.
How far down stream does this impact the coastal current? don't know that either-- but idealized runs with an estuary undergoing spring neap cycles should be assess to answer this. In the real world, spring neap modulation at the plume's source is likely to be eventually obscured by winds.
Bob Chant
The solution to the classic lock exchange (or internal dam burst) problem is that two fluids of differing density (rho1 in x<0, and rho2 in x>0) adjust by propagating internal waves in both directions.
In J.S. Turner's book "Buoyancy Effects in Fluids", he quotes the well known result that the minimum energy solution is that with the upper layer thickness, H, being half the total water depth, d. i.e. H/d = 1/2, and for this ratio the velocity of the advancing nose of the lighter fluid is given as U = 0.5*sqrt(g'd) where g' is the reduced gravity (g'= g*(rho1-rho2)/rho1 = (approx) 0.3 m/s^2 for freshwater over seawater.).
Turner goes on to note (in section 3.2.5) that the speed of advance of a nose of light water over a deep fluid i.e. H<<d is given by U = sqrt(2)*sqrt(g'H). He adds that the velocity changes almost linearly with the ratio H/d between these two limits (H/d ~0 and H/d = 1/2). After a little bit of algebra, I get that one can write this linear dependence as:
U = sqrt(2)(1-H/d)sqrt(g'H)
I've been making the case in this discussion that I expect a useful criterion for a sensible river injection strategy in ROMS is that the thickness of the upper fresh layer should be that which gives a speed of advance for the nose that is matched to the average current in the river. In this case, the layer of fresh water spreads at its own natural speed under the force of a baroclinic pressure gradient maintained by drawing the steady flow of the river in behind it.
For a river flow rate of Q entering into a grid cell of width L, this criterion would state that Q/L should match UH. My algebra skills aren't up to solving the resulting equation involving H^3/2 but I have plotted HU = sqrt(2)*H*(1-H/d)*sqrt(g'H) vs H for a variety of values of d. See the plot at http://marine.rutgers.edu/~wilkin/UHvsH.gif. If you know your river transport Q and cell width L then you can start on the y-axis at Q/L and go over to the curve, and down to see the value of H that fits Turner's relationship. The total water depth, d, doesn't make a heck of a lot of difference.
For the Hudson River at Q = 600 m^3/s entering into a cell of width L = 1 km, Q/L=0.6 m^2/s. For a water depth in the cell of 20 m, I get H = 1 m as the upper layer depth that would correspond to how the plume would likely look after undergoing the dynamical adjustment to the lock exchange process. This is presumably a choice that would not drastically violate the hydrostatic approximation of the model. Note that H = 1 m is one level of my 20-level model.
I'm attracted to this argument because it's based on physics.
I do have to add that I found it somewhat unsettling to read Rob's comment:
Remember the words of Nobel laureate Ernest Rutherford: "All science is either physics or stamp collecting."
John.
In J.S. Turner's book "Buoyancy Effects in Fluids", he quotes the well known result that the minimum energy solution is that with the upper layer thickness, H, being half the total water depth, d. i.e. H/d = 1/2, and for this ratio the velocity of the advancing nose of the lighter fluid is given as U = 0.5*sqrt(g'd) where g' is the reduced gravity (g'= g*(rho1-rho2)/rho1 = (approx) 0.3 m/s^2 for freshwater over seawater.).
Turner goes on to note (in section 3.2.5) that the speed of advance of a nose of light water over a deep fluid i.e. H<<d is given by U = sqrt(2)*sqrt(g'H). He adds that the velocity changes almost linearly with the ratio H/d between these two limits (H/d ~0 and H/d = 1/2). After a little bit of algebra, I get that one can write this linear dependence as:
U = sqrt(2)(1-H/d)sqrt(g'H)
I've been making the case in this discussion that I expect a useful criterion for a sensible river injection strategy in ROMS is that the thickness of the upper fresh layer should be that which gives a speed of advance for the nose that is matched to the average current in the river. In this case, the layer of fresh water spreads at its own natural speed under the force of a baroclinic pressure gradient maintained by drawing the steady flow of the river in behind it.
For a river flow rate of Q entering into a grid cell of width L, this criterion would state that Q/L should match UH. My algebra skills aren't up to solving the resulting equation involving H^3/2 but I have plotted HU = sqrt(2)*H*(1-H/d)*sqrt(g'H) vs H for a variety of values of d. See the plot at http://marine.rutgers.edu/~wilkin/UHvsH.gif. If you know your river transport Q and cell width L then you can start on the y-axis at Q/L and go over to the curve, and down to see the value of H that fits Turner's relationship. The total water depth, d, doesn't make a heck of a lot of difference.
For the Hudson River at Q = 600 m^3/s entering into a cell of width L = 1 km, Q/L=0.6 m^2/s. For a water depth in the cell of 20 m, I get H = 1 m as the upper layer depth that would correspond to how the plume would likely look after undergoing the dynamical adjustment to the lock exchange process. This is presumably a choice that would not drastically violate the hydrostatic approximation of the model. Note that H = 1 m is one level of my 20-level model.
I'm attracted to this argument because it's based on physics.
I do have to add that I found it somewhat unsettling to read Rob's comment:
To me, unphysical numbers (like salinity approaching 40 psu when we input FRESH water) set off alarm bells, and I don't care whether the model solution looks smooth or not. A criterion that amounts to "the model runs and the solution looks smooth so it must be OK" is never safe.I am less concerned about unphysical numbers as I am about
numerical instability when anything is resolved with only one grid cell.
Remember the words of Nobel laureate Ernest Rutherford: "All science is either physics or stamp collecting."
John.
- arango
- Site Admin
- Posts: 1361
- Joined: Wed Feb 26, 2003 4:41 pm
- Location: DMCS, Rutgers University
- Contact:
Well, I like scientific discourse. It is very clear that we need to test all these ideas. I am thinking that we need to set-up an idealized river runoff test problem in which we can explore all these issues.
Perhaps, we can have a composed grid application. In this way we can test also the nesting/composed grid capabilities of ROMS. What about a grid configuration like:
We can run this application with both grids (1,2) or single grid (1). If single grid, the river is modeled as a boundary condition. If composed grid, we can study the impact of the full estuary. In a composed grid, we can have different horizontal and vertical resolutions for grids 1 and 2. This application can also be configured with a single grid with land/sea masking, like the RIVERPLUME test. The advantage of a composed grid is that we can play also with the angle of grid 2 with respect grid 1. Then, we can compare the solutions. We can force with tides at the open boundaries of grid 1 and explore Bob Chant questions.
One thing that we need note here is that in river runoff we have two different processes: transport (mass flow) and salt wedge dynamics. This salt wedge dynamics is ignored when the rivers are treated as boundary conditions. There is a gravitational adjustment going here. Recall the classical gravitational adjustment between to fluids of different density. There is a lot of waves (gravity and internal waves) that can be generated at the sharp density front.
Hernan
arango@imcs.rutgers.edu
Perhaps, we can have a composed grid application. In this way we can test also the nesting/composed grid capabilities of ROMS. What about a grid configuration like:
Code: Select all
___
| |
| |
| |
| 2 |
| |
______|___|______
| |
| |
| 1 |
| |
|_________________|
One thing that we need note here is that in river runoff we have two different processes: transport (mass flow) and salt wedge dynamics. This salt wedge dynamics is ignored when the rivers are treated as boundary conditions. There is a gravitational adjustment going here. Recall the classical gravitational adjustment between to fluids of different density. There is a lot of waves (gravity and internal waves) that can be generated at the sharp density front.
Hernan
arango@imcs.rutgers.edu
Last edited by arango on Mon Apr 05, 2004 8:54 pm, edited 4 times in total.
Speaking from one who has spent many years studying river plumes, wanted to add the following...
I've only briefly skimmed this topic at Hernan's request, but I would have to say I completely agree with Hetland:
- Unless the surface layer (and I mean literally the layer thickness set by the model) is indeed where water is discharged in nature, there will likely be a gross error in dynamics both in transport and mixing of the buoyant fluid (more buoyant water exposed to surface effects, much higher inertia than in nature, very different Froude and Rossby number conditions).
- Rob puts it best in his first post:
Related to Enrique's suggestion, I'm sure one can clearly show that a discharge over the river mouth depth gives significantly different results than one discharge in only the surface layer (unless the model vertical resolution is such that the "top level/layer" approaches the depth of the river mouth. Perhaps, Rob could do two simple runs and put this issue to rest.
Derek Fong
I've only briefly skimmed this topic at Hernan's request, but I would have to say I completely agree with Hetland:
- Unless the surface layer (and I mean literally the layer thickness set by the model) is indeed where water is discharged in nature, there will likely be a gross error in dynamics both in transport and mixing of the buoyant fluid (more buoyant water exposed to surface effects, much higher inertia than in nature, very different Froude and Rossby number conditions).
- Rob puts it best in his first post:
Certainly, one can imagine cases where the surface layer is close enough to what the real river mouth discharge looks like. But in my experience, that would be more of an exception than the norm.Rivers are a boundary condition within the model. As such, they should be formed just like any other boundary condition -- in a way that is numerically sound and physically consistent.
Related to Enrique's suggestion, I'm sure one can clearly show that a discharge over the river mouth depth gives significantly different results than one discharge in only the surface layer (unless the model vertical resolution is such that the "top level/layer" approaches the depth of the river mouth. Perhaps, Rob could do two simple runs and put this issue to rest.
Derek Fong
If you read the above discussion carefully I think you'll see we are all largely in agreement. The apparent controversy stems from not acknowledging the different applications that people are addressing.
Derek Fong notes:
Something we do know for sure is that the lock-exchange test problem is a very challenging test of the numerics and severly prone to over/undershoot of the advection, and I believe accumulation of errors associated with this leads to the over/undershooting many of us have witnessed when running river source applications.
Derek adds that:
But there are a large number of us who are not so interested in the dynamics of an individual river inflow in detail as we are of capturing the broader scale impacts of riverine buoyancy inputs on shelf-wide circulation. We are therefore studying so-called ROFIs (regions of freshwater influence) with less regard for the details of how this buoyancy actually flows through small estuaries and bays to reach the shelf.
The shelf-wide modeling we are doing at Rutgers is in this category. We have a 9-km resolution model of the North American East Coast from Newfoundland to Cubu with 40 rivers. The 39th largest river is the Saco with an average volume transport of only 95 m^3/s (Our 40th river is the TOMS River in NJ which we just had to include even though there are other larger rivers ). The Saco River flow is therefore a trickle compared to the 9-km cell into which it must flow. If the inflows are done carelessly, the undershooting of the salinity in small underresolved estuaries is extreme, upsetting the equation of state and our biogeochemistry model.
Similarly, Enrique Curchister's US West Coast model has has relatively deep water at the coast necessitated by the steep shelf and stability and accuracy issues related to bathymetric steepness.
The query that sparked this discussion came from our Adriatic modeling colleagues, which is another major body of water with a ROFI shelf.
Rob states:
I'm only saying I put the rivers into the topmost layer because my math comes out to 1-layer for my application. Everyone has to reason this through for themselves for their rivers and grids.
I would like to second Derek's nomination...
John.
Derek Fong notes:
This is true. The discussion I offered on the lock exchange solution introduces some rigor into how to estimate what the inflow layer thickness would be to "pre-adjust" the inflow to the appropriate Froude number. A layer thinner than the H value I propose would presumably introduce a flow with excessive inertia, advancing faster than the internal wave speed. A thicker layer is subject to lock-exchange adjustment by generating a strong internal wave response.Unless the surface layer (and I mean literally the layer thickness set by the model) is indeed where water is discharged in nature, there will likely be a gross error in dynamics both in transport and mixing of the buoyant fluid (more buoyant water exposed to surface effects, much higher inertia than in nature, very different Froude and Rossby number conditions).
Something we do know for sure is that the lock-exchange test problem is a very challenging test of the numerics and severly prone to over/undershoot of the advection, and I believe accumulation of errors associated with this leads to the over/undershooting many of us have witnessed when running river source applications.
Derek adds that:
ROMS users like Derek and Rob Hetland who have studied river outflow and plume kinematics in detail have naturally configured their applications to resolve the processes they are interested in. They have presumably pushed their applications as close as possible to the limit of very shallow water depths and small grid sizes at the coast where inflows occur. Then the volume of the cell into which the river enters is modest compared to the flow rate times some characteristic time scale of interest, so that ROMS is able to perform the dynamical adjustment itself over several grid cells, effectively partially resolving the lock-exchange adjustment rather than parameterizing it.Certainly, one can imagine cases where the surface layer is close enough to what the real river mouth discharge looks like. But in my experience, that would be more of an exception than the norm.
But there are a large number of us who are not so interested in the dynamics of an individual river inflow in detail as we are of capturing the broader scale impacts of riverine buoyancy inputs on shelf-wide circulation. We are therefore studying so-called ROFIs (regions of freshwater influence) with less regard for the details of how this buoyancy actually flows through small estuaries and bays to reach the shelf.
The shelf-wide modeling we are doing at Rutgers is in this category. We have a 9-km resolution model of the North American East Coast from Newfoundland to Cubu with 40 rivers. The 39th largest river is the Saco with an average volume transport of only 95 m^3/s (Our 40th river is the TOMS River in NJ which we just had to include even though there are other larger rivers ). The Saco River flow is therefore a trickle compared to the 9-km cell into which it must flow. If the inflows are done carelessly, the undershooting of the salinity in small underresolved estuaries is extreme, upsetting the equation of state and our biogeochemistry model.
Similarly, Enrique Curchister's US West Coast model has has relatively deep water at the coast necessitated by the steep shelf and stability and accuracy issues related to bathymetric steepness.
The query that sparked this discussion came from our Adriatic modeling colleagues, which is another major body of water with a ROFI shelf.
Rob states:
The formulation of the rivers boundary condition I propose is motivated by physics, not a guess (well maybe a bit of a guess). The numerical treatment in ROMS imposes the river sources by augmenting the fluxes on the cell faces. Subsequent divergence of these fluxes adds to the tendency terms for the depth-averaged momentum and tracers. How the divergence in the source cell is subsequently mixed and transported away depends on the resolution and other mixing parameterizations used in the model. I assume that even with unpleasant over/undershooting of the salinity that the total divergence is being done correctly (first moment conservation) and the total buoyancy input is correct. Arguably, some cranked up mixing would obscure the problem and those of us who care only about the total buoyancy input in a ROFI could live in blissful ignorance of the problem. However, the potential for adverse impact on the (nonlinear) biogeochemistry means we need to give this a little more thought.Rivers are a boundary condition within the model. As such, they should be formed just like any other boundary condition -- in a way that is numerically sound and physically consistent.
I'm only saying I put the rivers into the topmost layer because my math comes out to 1-layer for my application. Everyone has to reason this through for themselves for their rivers and grids.
I would like to second Derek's nomination...
How about it Rob?Perhaps, Rob could do two simple runs and put this issue to rest.
John.
First of all, I was unclear when I said unphysical numbers. I meant local values of the outflow speed and Froude number not matching measured values but otherwise physically consistent, rather than numerical dispersion (i.e., the very high salinity). These slightly elevated values of salinity near the bottom are impossible to avoid completely in numerical advection schemes without flux correction. However, we can minimize the effects so that the problem essentially goes away. I measure this effect in terms of spurious fresh water formed. You can calculate the amount of (spurious) fresh water required to eliminate the spurious dense water formed (see attached paper). I should note that, even when you have a single grid point with very high salinity, it may not affect the solution much, if the salinity remains constant in time and does not begin to flow out into the rest of the domain.
The part missing from this discussion is what are the dynamical scales of the plume? We can only hope to simulate processes that we can resolve numerically (hence, my insistence that the outflow must be put over more than one point). This usually means 3+ grid points. The scales of the plume change, depending on which processes are locally important. This means as water is moved from the river to the ocean, there are many different dynamical scales, depending on the local physics. In my opinion, local physics is represented by the dominant entrainment mechanism: either advective shear mixing (the near-field with scales of 1-10 km), or wind mixing (the far-field, with scales of 10-100 km).
In your Hudson plume model, you are not resolving the near-field. There is probably no point where the outflow becomes supercritical near the estuary outflow (in fact -- no estuary). This model is obviously designed to look at the far-field region of the plume, and that's fine.
If you were resolving the near-field, the shelf scales would suffer, or not even be included in the domain. My results (in the paper below), suggest that the far-field may not even be very sensitive to the details of the near-field, so that this model may still be accurate, despite missing the near-field details.
That said, any model that does not resolve sub-gridscale processes should ensure that these small scales (grid-scale noise in the larger model) do not contaminate the larger scale (resolved) part of the simulation. There no good way to do this yet. I would suggest we need to proceed by first understanding the near-field physics, so that we can think of ways to parameterize those processes when they are unresolved. The same holds true for estuarine dynamics. Bob's comments about the residence time of fresh water in the estuary are dead on. Not only don't we deal with changes in estuarine mixing, we don't even account for the changes in fresh water outflow due to steady mixing in the estuary. (The estuary will tend to smear out a pulse of fresh water, even when mixing is constant. I've attached a paper that talks about that as well).
Attached is a paper on wind-forced river plumes that I recently submitted to JPO that deals with the different dynamical scales of the plume:
http://www.myroms.org/links/wind_plume.pdf
I hope some of you will find it helpful; comments are very welcome.
Also attached is a paper on numerical simulations of Hansen and Rattray style estuarine dynamics:
http://www.myroms.org/links/estuary_scales.pdf
I think this is relevant because this is the sort of estuary that we will be able to include in our larger scale shelf simulations.
Rob
The part missing from this discussion is what are the dynamical scales of the plume? We can only hope to simulate processes that we can resolve numerically (hence, my insistence that the outflow must be put over more than one point). This usually means 3+ grid points. The scales of the plume change, depending on which processes are locally important. This means as water is moved from the river to the ocean, there are many different dynamical scales, depending on the local physics. In my opinion, local physics is represented by the dominant entrainment mechanism: either advective shear mixing (the near-field with scales of 1-10 km), or wind mixing (the far-field, with scales of 10-100 km).
In your Hudson plume model, you are not resolving the near-field. There is probably no point where the outflow becomes supercritical near the estuary outflow (in fact -- no estuary). This model is obviously designed to look at the far-field region of the plume, and that's fine.
If you were resolving the near-field, the shelf scales would suffer, or not even be included in the domain. My results (in the paper below), suggest that the far-field may not even be very sensitive to the details of the near-field, so that this model may still be accurate, despite missing the near-field details.
That said, any model that does not resolve sub-gridscale processes should ensure that these small scales (grid-scale noise in the larger model) do not contaminate the larger scale (resolved) part of the simulation. There no good way to do this yet. I would suggest we need to proceed by first understanding the near-field physics, so that we can think of ways to parameterize those processes when they are unresolved. The same holds true for estuarine dynamics. Bob's comments about the residence time of fresh water in the estuary are dead on. Not only don't we deal with changes in estuarine mixing, we don't even account for the changes in fresh water outflow due to steady mixing in the estuary. (The estuary will tend to smear out a pulse of fresh water, even when mixing is constant. I've attached a paper that talks about that as well).
Attached is a paper on wind-forced river plumes that I recently submitted to JPO that deals with the different dynamical scales of the plume:
http://www.myroms.org/links/wind_plume.pdf
I hope some of you will find it helpful; comments are very welcome.
Also attached is a paper on numerical simulations of Hansen and Rattray style estuarine dynamics:
http://www.myroms.org/links/estuary_scales.pdf
I think this is relevant because this is the sort of estuary that we will be able to include in our larger scale shelf simulations.
Rob
John - my apologies for not reading the whole thread thoroughly before commenting. You're right. For large scale problems, surface layer discharge may be all that one can resolve within the grid.
Rob - if you can take the time to elucidate the effects, etc., I think many would find the results/comparison insightful.
I guess it all boils down to wishing we could all be doing DNS of the ocean, but we can't. So compromises, judiciously, must be made whether in grid resolution, boundary forcing, or size of domain.
Derek
Rob - if you can take the time to elucidate the effects, etc., I think many would find the results/comparison insightful.
I guess it all boils down to wishing we could all be doing DNS of the ocean, but we can't. So compromises, judiciously, must be made whether in grid resolution, boundary forcing, or size of domain.
Derek
Estuarine Boundary Condition Notes
I would like to amplify some of Bob Chant's thoughts. Estuaries take a little bit of river flow, pull in a LOT of ocean water, mix the two, and spit out the mixture, which is usually only slightly fresher than oceanic. Two points then are worth noticing when formulating a "river" BC: (i) the outflow is usually only a little bit fresh and it has a lot more volume flux than the river alone (typically a factor of 10), and (ii) there is an inflow in the lower half of the estuary mouth that is about the same volume flux. In my numerical experiments (like the ROMS riverplume setup, but with a much better-resolved estuary) this inflow gets pulled from the ocean along the right-hand shore, almost like a reversed riverplume coastal current. This was predicted by Beardsley and Hart, 1978, JGR, 83, C2, 873-. Also, as Bob Chant mentions, these properties have order-one changes over time, as tidal mixing varies.
Parker MacCready
Parker MacCready
Hi
I decided to tag this on to the rivers discussion; I have a question with a slight change of emphasis. I am just adding a river to what I hope will be a semi-realistic model of the Gulf of Cadiz. The only river data I've managed to obtain come from a dam that is ~150 km inland (Alcala del Rio station on the Guadalquivir River). This leaves me with a problem as to how to ensure that the measured discharge enters the Gulf with the correct time lag. Clearly this would be done with knowledge of the flowrate, but I am interested in knowing if other people have had experience of similar situations, and how they have dealt with them? Any rules of thumb?
Thanks
Evan Mason
I decided to tag this on to the rivers discussion; I have a question with a slight change of emphasis. I am just adding a river to what I hope will be a semi-realistic model of the Gulf of Cadiz. The only river data I've managed to obtain come from a dam that is ~150 km inland (Alcala del Rio station on the Guadalquivir River). This leaves me with a problem as to how to ensure that the measured discharge enters the Gulf with the correct time lag. Clearly this would be done with knowledge of the flowrate, but I am interested in knowing if other people have had experience of similar situations, and how they have dealt with them? Any rules of thumb?
Thanks
Evan Mason
Re: Treatment of rivers
Regarding the internal waves generation problem mentioned by John Wilkin in an early post on this thread... how difficult would it be to add an option, and modify the code, to set the river T and S to the ambient ocean values at each vertical level? Would that solve the problem? Seems like this would be useful for adding neutrally buoyant dyes (anywhere in the domain, not just at an estuary), for example.
From my reading of the documentation and forum, using "river_dye_XX" is the only way to add a dye-type (as opposed to a discrete "float") passive tracer to ROMS (if running parallel). (Well, aside from running an NPZ model and tracking phytoplankton or something with the loss/gain turned off...which seems pretty clumsy.) But at present you can't add a river_dye without specifying T, S, and the direction and mass transport of the influx. Or am I wrong?
From my reading of the documentation and forum, using "river_dye_XX" is the only way to add a dye-type (as opposed to a discrete "float") passive tracer to ROMS (if running parallel). (Well, aside from running an NPZ model and tracking phytoplankton or something with the loss/gain turned off...which seems pretty clumsy.) But at present you can't add a river_dye without specifying T, S, and the direction and mass transport of the influx. Or am I wrong?
Re: Treatment of rivers
If you have inflowing water via a point source, you need to specify all the tracers to avoid the unstable downwind advection scheme. You can specify a climatology for the incoming T and S.
If on the other hand you want your dyes to change as they come through a passage, it is possible to code that up. Someone at Rutgers passed me the code for a TRC_PSOURCE.
If on the other hand you want your dyes to change as they come through a passage, it is possible to code that up. Someone at Rutgers passed me the code for a TRC_PSOURCE.
Re: Treatment of rivers
Kate,
thanks for your reply. Now that I have actually implemented a source in them model interior, I understand a bit better how it works, and will put some hints on a more active and relevant thread. But as to what I was asking about, it turns out this was implemented with 3.5 using LtracerSRC (for example LtracerSRC = F F T forces the incoming temperature and salinity to match ambient values). (When I noticed this feature, documented in a ROMS wiki, I imagined Hernan rolling his eyes at my ignorance.)
Following your email, I looked at a few of the tracer advection schemes in ROMS and didn't see any that explicitly advertised itself to be "downwind" (maybe you meant a momentum advection scheme?), and in any case if you use TS_MPDATA (which governs both horizontal and vertical advection) (and TS_MPDATA seems to be the preferred choice of those that know what they are talking about on this forum) ...then you are using an upwind scheme. According to a wiki on upwind schemes, first order upwind schemes tend to be subject to "severe problems with numerical diffusion in the presence of a strong gradient" - so it might be advisable to ramp up the transport or concentration at the start.
...that notwithstanding, no doubt you are right about having to provide values for that ALL of the tracers (e.g. temperature and salt) even if just assigning them zeros.
Thanks again and I hope you are enjoying the beautiful Fairbanks May weather.
thanks for your reply. Now that I have actually implemented a source in them model interior, I understand a bit better how it works, and will put some hints on a more active and relevant thread. But as to what I was asking about, it turns out this was implemented with 3.5 using LtracerSRC (for example LtracerSRC = F F T forces the incoming temperature and salinity to match ambient values). (When I noticed this feature, documented in a ROMS wiki, I imagined Hernan rolling his eyes at my ignorance.)
Following your email, I looked at a few of the tracer advection schemes in ROMS and didn't see any that explicitly advertised itself to be "downwind" (maybe you meant a momentum advection scheme?), and in any case if you use TS_MPDATA (which governs both horizontal and vertical advection) (and TS_MPDATA seems to be the preferred choice of those that know what they are talking about on this forum) ...then you are using an upwind scheme. According to a wiki on upwind schemes, first order upwind schemes tend to be subject to "severe problems with numerical diffusion in the presence of a strong gradient" - so it might be advisable to ramp up the transport or concentration at the start.
...that notwithstanding, no doubt you are right about having to provide values for that ALL of the tracers (e.g. temperature and salt) even if just assigning them zeros.
Thanks again and I hope you are enjoying the beautiful Fairbanks May weather.
Re: Treatment of rivers
I wasn't kidding when I said "unstable downwind scheme". I have had a source in which I asked it to use the downwind value for the tracers and I watched it go bad in the debugger. Then someone who knows their numerics better than I do pointed out that sure, downwind schemes are known to be unstable. Use them at your peril!
Ah, Fairbanks in May! The leaves burst out yesterday, though they are still tiny. I just love the light. Why is it that I have two trips planned?
Ah, Fairbanks in May! The leaves burst out yesterday, though they are still tiny. I just love the light. Why is it that I have two trips planned?