Calculating sediment flux across transect

[mbiddle]

Hey All,

I'm trying to compute the total amount of suspended sediment that passes through a transect in my domain over the simulation time series. Originally, I was using the velocity data compiled with the grid shape and suspended sediment concentration to compute the fluxes. That seemed overly complicated. So I dug back into the manual and saw that sediment.in has a flag for tracer volume flux (mud and sand), how handy! :

Code: Select all

Aout(MHUTav) == T       ! Huon_mud_01, ...        tracer volume flux, <Huon*t>
Aout(MHVTav) == T       ! Hvom_mud_01, ...        tracer volume flux, <Hvom*t>

So, I set those to true and re-ran the simulation. Now that I have the flux data, I noticed they are on different grids (u and v grid) for the flux across the faces of the cells. Which conceptually makes sense.

However, the grid shapes are not the same size. Below is the shape of the flux variables for 2497 time-steps and 5 sigma layers (my grid is 100x100):

Code: Select all

f.variables['Huon_mud_01'].shape
(2497, 5, 100, 99)
f.variables['Hvom_mud_01'].shape
(2497, 5, 99, 100)

Initially I thought I could translate the u and v grids to something consistent, so I went digging around and found this function in pyroms. Which can translate u or v grid to psi grid. Here's my snippet for the translation:

Code: Select all

Huon_sand_01 = 0.5 * (Huon_sand_01[:,:,1:,:] + Huon_sand_01[:,:,:-1,:]) # u ==> psi
Hvom_sand_01 = 0.5 * (Hvom_sand_01[:,:,:,1:] + Hvom_sand_01[:,:,:,:-1]) # v ==> psi

The response has agreeable dimensions between the u and v values, so I can start performing the calculations to get the along and cross channel flux. But, from my understanding the psi grid is the node, not the center (rho) of the cell. So, I don't think this is the translation I need to perform.

Does anyone have guidance (or examples in python, matlab, or anything else?) on how to compute the flux across a transect? I think I've overthought this one...

Thanks in advance!

[kate]

Are you aware of how the whole grid looks? See for example the figure here.
To get something in the centers, try something like:

Code: Select all

Huon_sand_01 = 0.5 * (Huon_sand_01[:,:,1:-1,1:] + Huon_sand_01[:,:,1:-1,:-1]) # u ==> rho
Hvom_sand_01 = 0.5 * (Hvom_sand_01[:,:,1:,1:-1] + Hvom_sand_01[:,:,:-1,1:-1]) # v ==> rho

The result should be 98x98, the size of your interior computational rho grid (no image points).

[mbiddle]

Thanks, Kate! That figure was very helpful.

[wilkin]

The Huon_* quantities are defined on cell faces, i.e. respective u or v faces, because this is the way ROMS does fluxes internally. The divergence in "x" of fluxes on two neighboring "u" faces gives the cell centered "rho" (or "sand") point flux. So following a wiggly path of u and v faces will give you the most precise net flux through a transect.

[mbiddle]

Okay, I adjusted my approach. I'm trying to see if there is at least a balance in the water flux over the grid. So I started investigating the Huon and Hvom variables, ignoring the sediment ones for now.

Instead of translating the grids, I think it's okay to use the existing face values on their respective grids. If I understand correctly huon(x,y) is the flux across the u face at grid index (x,y) and hvom(x,y) is the flux across v face at the same grid index (x,y).

If I have the indexes for the transect from which I'd like to extract the flux data, I can simply extract those values by slicing the huon and hvom data by the (x,y) indexes.

To calculate the total flux across a transect, I calculate the sum of the Huon flux and the sum of the Hvom flux (horizontally and vertically) for each point in time. Now I have one Huon and one Hvom at each point in time for each transect. I can then convert from North-South fluxes to across and downstream fluxes (since the flux isn't north-south). This gives me a time-series of the across and downstream flux at each transect.

I would assume, with my model configuration of a major river input source at the north-west boundary, that there would be a balance of fluxes across each transect. The water coming in at the river would be relatively equal to the water exiting the system (at the southwest boundary).

I've attached figures depicting the raw flux time-series, the determination of the angle of rotation, and the rotated flux values for each transect. I've also included a map of where the transects are in relation to the model's mask_rho.

Tb - northern boundary transect
T1 - middle transect
T2 - southern transect

The problem I'm having is the water fluxes aren't balancing between all three transects:
Tb max flux = 23509.4 m3/s
Tb integrated flux = 5.48596e+06 m3/s

T1 max flux = 23586.4 m3/s
T1 integrated flux = 5.38838e+06 m3/s

T2 max flux = 18894 m3/s
T2 integrated flux = 4.0243e+06 m3/s

Tb and T1 are close enough in both the maximum flux and the time integrated flux. However, T2 is significantly lower. Basically, this is telling me that water is coming into the region, but not exiting. Where is it all going?

The east boundary has a land mask blocking any water flow, this is similar for the north-east boundary. I do have a a varying free surface for the southern boundary to act as tidal input. But, I don't think that would cause this reduction in the fluxes.

I've tried to do a stair-step approach for transect T1 and set it up to be perpendicular to the flow. But this resulted in a really low total flux across that transect (max flux = 12418.6 m3/s), after performing the same calculations.

Could you please elaborate on

So following a wiggly path of u and v faces will give you the most precise net flux through a transect.

Any help would be greatly appreciated!