gdal translate - ETOPO1 upsampling issue

I am using the OSGeo4W Shell with the ETOPO1, 1 arc-minute global bathymetric data (etopo1_ice_g_f4.flt).

I am cutting 256x256 pixel slippy map tiles, and am having a problem with getting smooth-appearing PNG files out of the GDAL process. For example, the tile at zoom level 10, x = 306, y = 389 represents a 24 arc-minute span, from 39.24 N to 39.64 N, and 72.42 W to 72.02 W.

So, to produce the GeoTIFF I apply:

gdal_translate -of GTiff -projwin -72.42 39.64 -72.02 39.24
-outsize 256 256 -r bilinear
F:\test\etopo1_ice_g_f4.flt F:\test\z10x306y389.tif

Then, to get the color relief tile, I apply:

gdaldem color-relief 

F:\test\z10x306y389.tif F:\test\colors.txt F:\test\z10x306y389.png 
-of PNG

The problem is that the resulting PNG file does not appear to have been properly interpolated. Instead, it is pixelated:

I have asked a similar question, thinking that the problem was because I was using a .NET wrapper for the GDAL libraries, but here I am using OSGeo4W Shell, and the results are the same: a pixelated image. So, this question is more about upsampling ETOPO1 data.

The pixelation remains the same regardless of the type of resampling method I choose (e.g., bilinear, cubic, lanczos, etc.), so I wonder if this is an issue with the ETOPO1 data source?

The intermediary TIFF file (z10x306y389.tif in my example) opens as "all black" in various image editors (e.g, Paint).

And the gdaldem hillshade option results in a simply terrible graphic:

Has anyone else had issues when upsampling ETOPO1 elevation data? Is there a way to check the "correctness" of the steps I am using, or the integrity of the intermediary GeoTIFF file?