Your approach essentially has two major problems:

• The USGS 1 arc second files contain bad quality void fill in some tiles that would need to be detected and removed for good quality results.
• The void filled 3 arc seconds data from the USGS/NGA vary strongly in quality of the void fill, in Europe they mostly just interpolate (slightly better than your gdal_fillnodata attempt but essentially similar).

I agree with this; they’re quite evident when you see the hillshading layer.

You can see that for example in your final image where the valley at the bottom appears dammed near Mollières.

If by dammed you mean that smooth looking slope, yes, of course, I’m using less accurate data to fill that void. What I’m saying that this is the better I get from what I have.

Here for comparison a rendering of my current basic processing in that area

I would like to know which data are you using. I downloaded deFerrantis 1 arc for that region (http://www.viewfinderpanoramas.org/dem1/N44E007.zip) and indeed they look good in this region, but if you check the plains to the NE, you’ll see not nice terracing, like the one I showed in my post. There are also several vertical and horizontal ‘cuts’. Just run gdaldem hillshade on it to see what I mean.

I prefer to balance a bit of accuracy for time saved (again, I’m not a professional cartographer like you, this is mostly for having fun… and useful offline maps in my phone) and better looking plains.

with OSM lakes flattened.

Well I don’t need this, as I’m drawing the lakes on top of the terrain anyways. See:

http://grulicueva.homenet.org/~mdione/Elevation/#13/44.1810/7.1582

(This is not the definite map; in that one, slope shades are too dark and confusing).