First off, excellent job on creating a minimal working example. That really helps me understand where you're at and what you're trying to do. I was able to create this plot using your code:
We'll tackle your Point #2 first: When you try to plot color-scaled circles in the same map as color-scaled contour lines, Matlab gets confused. That's not your fault. That's Matlab's fault. The ability to plot two colormaps on the same plot is something users have been requesting for decades, but for some reason Matlab is still refusing to allow it. There are a few workarounds, but they're sort of clunky and probably not worth getting into here.
For your case, I recommend that you either use color for surface elevation, or you use color for the scattered data. Here's how you could use color for the scattered data and make all of the contour lines gray:
Lat = linspace(-76,-75,30)';
Long= linspace(120, 150, 30)';
scatterps(Lat, Long,50, Sodium,'filled')
graticuleps(-88:2:-76,-150:30:180)
h = bedmachine('surface','contour',0:100:5000,'k');
caxis([min(Sodium(:)) max(Sodium(:))])
If you want to label those contours, you'll have to take a manual approach. Instead of simply calling bedmachine to plot the data, you'll have to load the bedmachine data using the x and y limits of the map, then contour and label it, like this:
[sfz,x,y] = bedmachine_data('surface',xlim,ylim);
[C,h] = contour(x,y,sfz,0:100:5000,'k');
caxis([min(Sodium(:)) max(Sodium(:))])
clabel(C,h,'labelspacing',1000,'fontsize',7,'color',0.6*[1 1 1])
You could use other data like a background image for better context:
modismoaps('contrast','low')
Or if you want to use the size of circles to represent data values, here's how you can do it (caveat below).
Lat = linspace(-76,-75,30)';
Long= linspace(120, 150, 30)';
circleps(Lat, Long,Sodium,'facecolor',[.01 .26 .87],'facealpha',0.8)
graticuleps(-88:2:-76,-150:30:180)
bedmachine('surface','contour',0:100:5000);
[lat_legend,lon_legend] = ps2ll(xl(1) + diff(xl)/10,yl(2)-diff(yl)/15-(1:5)*diff(yl)/15);
circleps(lat_legend,lon_legend,5:10:45,'facecolor',[.01 .26 .87],'facealpha',0.8)
textps(lat_legend,lon_legend,{'5';'15';'25';'35';'45 ppb'},'horiz','center','vert','middle','color','w')
In the map above, the circles are sized based on 1 km radius per 1 ppb Sodium. And you could easily multiply the circle radii by some factor to make them visible on your map.
A word of caution about this approach though: On the face of it, it makes sense to scale the size of circles with the numeric value of your Sodium data, and this type of bubble chart is very common in geospatial data visualization, especially in popular media. However, Edward Tufte opened my eyes to the fact that size-scaled bubble charts are intrinsically inaccurate.
In the example above, I've scaled the circle raddi linearly with the Sodium data, but of course the area of a circle scales as the square of the radius. In this sense, the linear scaling of the radius exaggerates the size of the large data and minimizes the small data.
It might be better to scale the circle radii with the square root of the Sodium data, because then the amount of ink on the page devoted to each circle will scale linearly with the underlying data. However, the human brain does not perceive the size of circles as a linear function of radius or area. In truth, we tend to perceive circle sizes somewhere in between.
There's a feature in ArcGIS that acknowledges the fact that neither radius nor area is a good way to scale numeric data, and so in ArcGIS you can select an option to plot circle sizes based on some approximate human perception, but that introduces yet another problem, because then it becomes impossible to measure any part of the visualization and tie it directly to the underlying data.
Anywho, I try to stick with color scaling for numeric data, and single-color labeled contours and/or hillshade when plotting surface elevation or context. But do whatever makes most sense for your data.