My method in Terra3D was to generate random height values, and then make 3 smoothing passes with a little reshaping written into the algorithm. The reshaping causes the terrain under the water line to shift downwards a bit and everything above the water line shift up a bit. The effect though is a lot of hills and small lakes. And that may not be what you're looking for exactly.
However, I'm not convinced that the principal behind my method couldn't be used to get what you want. You may just have to write some more conditions into the reshaping/smoothing algorithm. For example, reducing the amount of smoothing on the terrain at higher elevations will create more of a rocky mountain appearance. And then writing a more involved smoothing (or gaussian) algorithm that stretches out further for the lower elevations could pull those lakes together to form natural rivers.
Here's the code for the terrain generation in Terra3D in case you're interested:
// GENERATE TERRAIN
for (i = 0; i < MAX; i++)
{
for (i2 = 0; i2 < MAX; i2++)
{
if (i<10 || i2<10 || i>MAX-10 || i2>MAX-10)
field[i][i2].y=0;
else
field[i][i2].y=(GLfloat(rand()%151)-75)/50+(field[i-1][i2-1].y+field[i-1][i2].y+field[i-1][i2+1].y+field[i-1][i2-2].y+field[i-1][i2+2].y)/5.05;
}
}
// SMOOTH/RESHAPE TERRAIN
for (int cnt = 0; cnt < 3; cnt++)
{
for (int t = 1; t < MAX-1; t++)
{
for (int t2 = 1; t2 < MAX-1; t2++)
{
field[t][t2].y = (field[t+1][t2].y+field[t][t2-1].y+field[t-1][t2].y+field[t][t2+1].y)/4;
if (cnt == 0)
{
if (field[t][t2].y < -1 && field[t][t2].y > -1-.5) field[t][t2].y -= .45, field[t][t2].y *= 2;
else if (field[t][t2].y > -1 && field[t][t2].y < -1+.5) field[t][t2].y += .5, field[t][t2].y /= 5;
}
}
}
}
It's sloppy code that I wrote about 10 years ago. Terra3D was a long learning process of experimentation and fumbling in the dark to produce the type of effects I was looking for. But maybe it'll help.