Melt rate vs u* and melt rate vs thermal driving

[claireyung]

Related to https://github.com/misomip/isomip-plus/issues/9

Scatter plots of these values could inform us which of the terms within the melt parameterisation are most important

[claireyung]

@dgwyther and Qin:

At the second hackathon we also revisited the idea of scatter plots of m/gamma_t vs uT and m/gamma_s vs uS . We decided these would be useful as a guide to the nonlinearity in the melt rate e.g. time averaging eddies and a verification tool. Would you be able to do these as well?

[dgwyther]

@claireyung what do you mean by uT and uS? Is this top layer velocity x temperature/salinity? Or maybe u* x temp or salinity?

[xylar]

@dgwyther, this was at least partly Ben's suggestion. I think it was u x T and u x S (friction velocity times either thermal or haline driving). A plot of m/Gamma_T vs. u x T would be expected to be a straight line for each model (see equation (24) in Asay-Davis et al., 2016). u x S is more complicated, see equation (26), because that is proportional to m * S_zd and S_zd isn't available in our output.

[dgwyther]

@xylar @claireyung Are the model GammaT settings stored in a netcdf file, or do they exist only in the PDF description for each model? edit: don't worry, I've just added them as a dictionary

[xylar]

I believe they are only in the pdf. They are also in @claireyung's table on Overleaf.

[dgwyther]

Just a note that there is a sizeable difference in the computed mean metrics, depending on the order in which masking and averaging is done:

melt = (data.meltRate*365*24*60*60).isel(nTime=slice(-13, None)).mean('nTime')
melt_av1= (melt.where(melt<1e36)).mean(dim='nx').mean(dim='ny')

produces a different value to:

melt = (data.meltRate*365*24*60*60)
melt_av = (melt.where(melt<1e36)).isel(nTime=slice(-13, None)).mean('nTime').mean(['nx','ny'])

I would have naievely thought they would produce the same fields, but they don't (at least in my testing). I would have thought that the order of the masking shouldn't matter, as in both cases (masking before or after the time mean), the masking will remove the bad cells. But possibly it's also behaviour of xr.where which I'm not clear about.

The order should probably be standardised for other plots.

[QinZhouLi]

I have pushed the changes under fork and branch: https://github.com/QinZhouLi/isomip-plus/blob/melt_vs_drivers/notebooks/melt_vs_drivers/melt_vs_drivers.ipynb.

Below are the figures of melt vs T and u_ for Ocean1_com:

and Ocean2_com:

Below are the figures of scaled melt (by constants and Gamma T) vs T times u_ for Ocean1_com

and Ocean2_com

According to the equation 24, the dots should align with the black dashed line (slope = 1). Surprisingly, most of the models are not.

[xylar]

@dgwyther, I think only the second one (taking the mean simultaneously in x and y) is correct and that is what we should use everywhere. If you take the mean first in x and then in y, the mean in y incorrectly assumes that each y value should get equal weight in the mean (unless there are no valid values at a given y value). Instead, each y value should have been given a different weight depending on how many valid values in x there were. So it makes sense to me that they given different results.

There's a simple test:

#!/usr/bin/env python
import numpy as np
import xarray as xr

ds = xr.Dataset()
ds['x'] = ('x', np.linspace(-2., 2., 101))
ds['y'] = ('y', np.linspace(-2., 2., 101))
ds['r'] = np.sqrt(ds.x**2 + ds.y**2)
ds['mask'] = ds.r < 1

ds['field'] = (1.0 - ds.r**2).where(ds.mask)

print(ds.field.mean(['x', 'y']).values)
print(ds.field.mean('x').mean('y').values)

The result is:

0.50573931
0.45691429

[dgwyther]

Ah hah! Yes that makes much more sense than the ordering of the masking. Thanks for figuring that out!

[xylar]

@QinZhouLi, that's interesting and troubling! We should try to figure out why most models aren't lining up right! Could you tell us which are lining up and which aren't? Among those lining up, I see:

• COCO
• MITgcm JPL
• POP2x

[QinZhouLi]

@xylar, both ROMS and FVCOM are also lining up.

[xylar]

Okay. My first guess would be that we have incorrect Gamma_T values for the other models, but that might be hard to prove.

[claireyung]

I have pushed the changes under fork and branch: https://github.com/QinZhouLi/isomip-plus/blob/melt_vs_drivers/notebooks/melt_vs_drivers/melt_vs_drivers.ipynb.

Thanks @dgwyther and @QinZhouLi ! Would you be able to make a pull request to the main repository too?

The scatter plots of melt vs T and u are interesting. Looks like both are noisily correlated to melt for Ocean1, less so for Ocean2, especially u which just collapses down to the minimum value of c_d^1/2u_tide. What happens if you plot for all points (and all time?) for a more local analysis rather than the whole cavity average - is it too many dots?

Okay. My first guess would be that we have incorrect Gamma_T values for the other models, but that might be hard to prove.

I compared the values you used, to what I have in the overleaf for COM experiments. I found these discrepancies:

GammaT={'COCO':2.5e-2,
'FVCOM':.2,
'MITgcm BAS':1.9e-2,
'MITgcm BAS CoupledV3':2.1e-2, 
'MITgcm JPL':3.25e-2,
'MOM6':4.5e-2, *-> should be 1.4e-1, but this is also the MOM6 COM experiment due to missing data so can ignore
'MOM6 SIGMA ZSTAR':8.7e-2 *-> I think this should be 4.5e-2 but need to check versions with Gustavo, 
'MPAS-Ocean':1.46e-2, *-> should be 1.94e-2
'NEMO-CNRS':4e-2, *-> should be 2.6e-2
'NEMO-UKESM1is':1.4e-1, *-> should be 4.5e-2
'POP2x':1.1e-1, *-> should be 1.1e-2
'ROMSUTAS':5e-2}

This might resolve the NEMO, MPAS and MOM6 ones but leaves the two MITgcm BAS ones not lining up with the 1:1 line. It also might change POP from kind of lining up to not really lining up... would you be able to run again with the updated values, @QinZhouLi, and confirm if that fixes things?

[QinZhouLi]

@claireyung , after updating those GammaT values, the figure of scaled melt (by constants and Gamma T) vs T times u_ for Ocean1_com is

and for Ocean2_com is

Now, there are only four models that don't line up: MITgcm BAS, MITgcm BAS CoupledV3, MOM6, and POP2X.

[xylar]

The value I have for POP2x in the PDF on my external hard drive is: Gamma_T = 0.0465085

The date at the top of my PDF is a little more recent than the one on my GoogleDrive so I think this might be correct.

Ocean1-4_COM_POP2x.pdf

@QinZhouLi, can you see if that value makes POP2x fall along the expected 1:1 line?

If so, we should replace the PDF in @claireyung's archive wh this one.

[QinZhouLi]

Updated figures with GammaT = 0.0465085 for POP2X for Ocean1_COM: and Ocean2_COM: It looks like the updated GammaT value is closer to the true value, but still not the same.

[xylar]

That's very odd. Thank you!

[xylar]

@QinZhouLi, the value in the PDF that I have on my Google Drive (rather than my external hard drive) is: Gamma_T = 0.1146 I believe that's what you tried originally and what looked good, right? Let's go back to that value.

@claireyung, I think the value of 1.1e-2 came from the TYP simulations, not the COM ones. The correct PDf is here. Could you replace the POP2x PDF with this one if they're not the same? Ocean1-4_COM_POP2x.pdf

[QinZhouLi]

@xylar, yes, with Gamma_T = 0.1146, POP2x looks much better, but still with some outliers, as Ocean1_COM: and Ocean2_COM:

[xylar]

@QinZhouLi, thanks! I think that's as good as we're going to get. Maybe the noise in the POP2x melt rates are associated with stronger nonlinearity or something. Not sure.

[claireyung]

Thanks @xylar and @QinZhouLi ! Indeed I had been accidentally reading the TYP transfer coefficient for POP2x rather than COM, so the PDFs are correct.